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ABSTRACT 


The subject of this investigation is to answer the question: 
Whether any significant increment to accuracy could be transferred from 
a super-control continental net (continental satellite net or super- 
transcontinental traverse) to the fundamental geodetic net (first-order 
triangulation). This objective was accomplished by evaluating the 
positional accuracy improvement for a triangulation station, which is 
near the middle of the investigated geodetic triangulation net, by 
using various station constraints over its geodetic position. 

This investigation on a 1858 kilometer long triangulation chain 
shows that the super-control net can provide a useful constraint to 
the investigated geodetic triangulation net, and thus can improve it 
only when the accuracy of super-control net is at least 1 part in 
500,000. 

The preliminary accuracy of super-transcontinental traverse is 
already better than this limiting accuracy of 1 part in 500,000; how- 
ever, the preliminary accuracy of continental satellite net is lower 
than this limiting accuracy of 1 part in 500,000. As such, continental 
satellite nets do not seem to provide any useful constraint, at least 
to this particular investigated triangulation chain. 
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1. INTRODUCTION 


Geodetic triangulation has been accepted as an accurate method of 

determining "precise" coordinates for the triangulation stations of relatively 

short chains. This well-accepted idea was also given in an article "How 

accurate is First-Order Triangulation?" [Simmons, 1950, pp. 53-561 with 

the following introductory words: 

The question is often asked, "How accurate is the position 
of a triangulation station, " or "To what accuracy are the 
distances between triangulation stations known?" These 
questions are difficult to answer, principally because 
first-order triangulation gives the optimum accuracy 
in the measurement of great distances and there is at 
present no super yardstick to which it can be compared. 

Two modern technological advancements, namely, satellites and 
electronic distance measuring (EDM) instruments, have questioned the 
first-order triangulation accuracy, especially if triangulation is extended 
to distances longer than 1000 km or more. In such extended triangu- 
lation systems systematic errors like lateral refraction, propagation of 
observational errors, residual polar motion effects on latitude, longitude 
and azimuth, etc. [Mueller, 1969, pp. 61, 86-87; Pellinen, 1970, pp. 34-35; 
Wolf, 1950, pp. 117], which cannot be eliminated, accumulate. Lately 
the question has been raised whether any significant increment to accuracy 
is "cascaded" from a 1:1 million 1000 km net through a 100 km net to 
local control over 10 km distances. 

The satellite triangulation and super-transcontinental traverse, being 
of the highest achievable accuracy of today, i.e., super-control net of 
"zeroth" order, constitute a modern geodetic super structure, within 
which the classical geodetic triangulation is supposed to provide a geodetic 
control densification. 

According to the classical geodetic concept a lower order system should be 
tied to a higher order system [Jordan/Eggert/Kneissl, 1958, Vol. IV. 1, p. 112]. 
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Statistically, this means that the variance-covariance of the higher order 
system, as a lower limit for accuracy, be at least compatible with the 
internal precision of the lower order system. For all practical reasons 
the accuracy of the higher order systems should be substantially better 
(by a factor of two to three) than the subordinated system, thus supplying 
a rigorous constraint in the reduction of the lower order system [Schmid, 
1969, p. 4]. 

The objective of this investigation is to answer the question: 

Whether any significant increment to accuracy could be transferred 
from a super-control net to the basic geodetic net (first-order triangu- 
lation). This objective was accomplished by evaluating the positional 
accuracy improvement for a triangulation station, which is near the middle 
of the investigated geodetic triangulation net, by using various station 
constraints over its geodetic position. 
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2. DATA AND ACCURACY ESTIMATES 
2.0 Data 

For the purpose of the present investigation, the triangulation of the 
western-half of the United States has been considered, as this is more 
accurate than that of the eastern-half of the United States [Simmons, 

1950, p. 54]. The investigation is done on the chain from Moses Lake, 
Washington to Chandler, Minnesota (Figure 2.0-1), as these two stations 
are also common on both the continental satellite net (CSN) and the 
super-transcontinental traverse (STT). The data used was supplied by 
the Triangulation Branch of Geodesy Division, and the Geodetic Research and 
Development Laboratory, both of the National Oceanic and Atmospheric 
Administration, Washington. 

2.01 Geodetic Triangulation Net. 

The observed data used is for the period 1897-1965. The approximate 
coordinates used are from 1965 adjustment. The details of Moses Lake - 
Chandler triangulation chain are as follows: 


Number of stations 


191 

Number of bases |~ Taped 


27 

LGeodimeter 


2 

Laplace stations 


13 

Observed directions 


919 

Distance between two stations 

Minimum 

273 m 


.Maximum 

190 km 

Total length of the chain 


I 858 km 
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The supplied data was in U. S. sign conventions, i.e., longitude 
positive westwards and azimuth reckoned clockwise from south [Bouchard 
and Moffitt, 1964, pp. 94, 315; Mitchell, 1948], which was converted 
internally within the program to the conventional sign conventions, 
i.e., longitude positive eastwards, azimuth reckoned clockwise from 
north [Grossmann, 1964, pp. 5-6; Mueller, 1969, pp. 15-19] . It is 
assumed that the necessary reductions have been applied to the observed data 
and the weight function P is "a priori" known to be a sufficient good accuracy . 

The configurations and specifications for triangulation net are dealt 
with in [Bomford, 1965; Bouchard and Moffitt, 1964; Gossett, 1950]. 

A typical configuration of U. S. Coast and Geodetic Survey triangulation 
chain is shown in Figure 2. 0-2. 

2.02 Super-Control Nets. 

Super-transcontinental traverse (STT) runs across the western-half 
and the eastern-half of the U. S.A. (Figure. 2. 0-3). Its specifications, con- 
figuration, reduction of data and instrumentation are dealt with by Meade 
Q967; 1969a; 1969b]. 

Continental satellite net (CSN) is, in general, planned in such a way 
so that the stations are around 1200 km apart and that these stations are 
evenly distributed over the entire area. CSN-stations are either the 
stations of first-order triangulation net or these are connected to them. 

Its specification and configuration are dealt with in [Deker, 1967; 

Mueller, 1964; Pellinen, 1970; Schmid, 1970; Shedlikh, 1970; Veis, I960]. 

The continental satellite net of North American Continent comprises of 
twenty stations which can be anchored in the three world net stations; 

Thule, Greenland,, Moses Lake, Washington,, and Beltsville, Maryland. 
Furthermore, planned is a tie to a fourth world net station - Shemya 
(Figure 2.0-4) [Schmid, 1970], 
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Triangulation Configuration. 
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Figure 2.0-4 Continental Satellite Net of North America 
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2. 1 Accuracy Estimates . 


2.10 Data . 

The following representative standard errors for observed data of Moses 
Lake - Chandler triangulation chain has been suggested ipMeade, >1)970]: 


Directions 0"4 


Azimuth 


0."8 


Base 


Taped 


L Geodimeter 


1 part in 500, 000 
1 ppm for distance > 15 km 
.1.5 cm for distances up to 15 km 


2.11 Networks . 

The mean of all section closures, which is the accuracy measure for the 
investigated geodetic triangulation net, is given as 1 part in 317,000 r Adams, 
19301. The standard position errors of the end stations of super-transcontinental 
traverse, which represent its accuracy measure, using actual data sets as given 
by different investigators differ too much from each other. The proportional 
error, which is the standard position error divided by the distance of the 
station from traverse-origin, is used for this investigation. The proportional 
errors of super-transcontinental traverse are given as follows: 1:740,000 over 
318 kilometer long traverse and 1:1,100,000 over 1270 kilometer long traverse 
TForeman, 1970]; 1:670,000 over 270 kilometer long traverse [Gergen, 1970] 
and 1:3,000,000 over 1858 kilometer long traverse TESSA., 19691. The pre- 
liminary accuracy (i. e. , proportional error) of continental satellite net, as 
obtained from the supplied data, corresponds to 1:385,000 for Chandler station. 
Because of this wide range in preliminary accuracy measures of these two super- 
control nets, investigations using the following accuracies (station constraints), 
are made: 1:300,000; 1:400,000; 1:500,000; 1:600,000; 1:700,000; 1:1 M; 

1:1. 5M; 1:3M. The use of these accuracy measures, which are within the limits 
of preliminary accuracies of the two super-control nets, will determine a limit on 
the accuracy requirement of the super-control net, which would be necessary to 
improve the geodetic triangulation net. 


- 9 - 





3. COMPUTATIONS AND RESULTS 

During the earlier period of this investigation considerable thought 
was given to the selection and use of such formulas and methods which 
would not only provide high accuracies, but also minimize or eliminate 
loss of accuracy in computations. This resulted in using Helmert - 
Rainsford - Sodano’s Iterative Solution for Inverse Problem (Appendix I), 
which are equally applicable for short and long lines, and Conjugate 
Gradient Method (Section 3.2) for the adjustment of the triangulation 
nets, where the original observation equation coefficient Matrix (A-Matrix) 
is used, thus avoiding direct formation of normal equations where certain pro- 
perties of the original A-Matrix are lost. To minimize the round-off errors, 
computations are done in double-precision with double-precision storage 
rMulIer-Merbach, 1970; Section 3. 28], 

3.0 Selection of Adjustment Method 

From the two basic adjustment methods, i.e., Method of Observation 
Equations and Method of Condition Equations, the former has been preferred 
for the present investigation due to reasons of simplicity and clarity. The 
reasoning of this preference has been dealt with in many publications 
("Ashkenazi, 1967, p. 167; Gotthardt, 1968, p. 180; Grossmann, 1961, 
p. 174; Helmert, I. Teil, 1880, p. 556; Jordan/Eggert/Kneissl, Bd. IV. 1, 
1958, p. 537; Wolf, 1968, p. 323; Wolfrum, 1969, p, ll. Due to the large 
size of the triangulation net under investigation and the availability of digital com 
puters, iterative methods were considered because (1) they are easier to program 
(2) they require less storage space as the coefficient matrix of a triangulation 
net is very sparse, (3) they use directly the original set of equations through- 
out the process and hence rounding-off errors do not accumulate from one 
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Iterative cycle to another. One important factor for deciding to use an itera- 
tive method is to know in advance whether the rate of convergence is rapid 
for the system. However, if the rate of convergence is found to be slow, which 
in general is the case in large systems, it can be accelerated bv a certain 
process T Ashkenazi, 1969, p. 34; Fox, 1965, p. 194; Hilger, et al. 1967, p. 11 
Ralston, 1965, p. 437], While searching for a suitable adjustment method this 
investigator came across the Conjugate Gradient Method (Cg - Method) f Ralston. 
1965; Schwarz, 1968 and 1970; Wolf, 1968; Zurmuhl, 19581. which has the 
following distinct advantages over other iterative methods, such as Gauss- 
Seidlr, Jacobi-, ^Relaxation- and. Gradient methods: 

1. Original A-Matrix is used, thus avoiding the formation of normal 
equations, where certain useful characteristics of A-Matrix. such as 
very small coefficients may be lost. 

2. Original A-Matrix, which has very few non-zero elements, is 
easily stored in comparatively much less computer space using an 
Index-Matrix. 

3. No "mesh-point numbering technique" TAshkenazi, 19671 to keep 

the band-width of the system a minimum is necessary. Thus stations 
can be added or taken out from the existing tri angulation system with- 
out caring for their numbering. 

4. It is a finite iterative process. Theoretically, the solution vector is 
obtained in a maximum of n-steps, n being the number of unknowns. 
Therefore, eigenvalues need not be calculated for determining the 
convergence. However, experimentation (Section 3. 4) shows that the 
solution vector is not obtained in n-steps, as the orthogonality between 
the residue-vectors is not maintained exactly (Section 3. 27). Con- 
sequently, the residue- vector r ( D ) after n-iterations is not zero 
(Section 3. 29). This deviation from zero depends upon the condition of 
the system; the poorer the condition, the larger will be the deviation. 

5. Even in case of a poorly conditioned system solution vector is obtain- 
able after more iterations. 
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6. In large nets the effect of round-off errors depends upon the elimination 
sequence in normal matrix [Korhonen, 1954]. In Cg- Method the 
elimination sequence plays no role since adjustment is simultaneous, 
i.e. , entire data as a whole is used. 

7. Each approximation x^ J ) to the solution vector is closer to the true 
solution x than the preceeding one. 

8. The ability to start anew at any point in the computation using the last 
x^ J ) as initial value so as to minimize the effects of round-off errors. 

3. 1 Mathematical Model (Method of Observation Equations) . 

Let Lj be the m independent observed quantities, v t the residuals to the 
observed quantities (obtained from the adjustment) and x,y, z, . . . the n unknown 
parameters to be determined. Each observation gives an observation equation, 
whose general from is 

Li + v t = f t (x, y, z, ...), (3.1-1) 

where i - 1,2,3. . . , m and f represents a linear or non-linear function. The 
method of least squares however demands that (1) f should be linear, i. e., 
a linear relationship between the observations and the unknowns and (2) the 
number of observations (m) should be greater than those of the unknowns (n) 
i. e.,m>n. In case of a non-linear function f this is linearized by using 
Taylor series about such good approximate values of the unknowns Xq, y 0 , z 0 , 
such that the second and higher order terms can be neglected. In this 
case, equation (3.1-1) can be written as 

Li + Vi = fi(Xo+dx, y 0 +dy, z 0 +dz, ... ) 
or 

Lj + Vj = fi(x 0 , y 0 , z 0 , ...)+ aidx + bidy + Cidz+ ... (3.1-2) 

i.e. v { = aidx + b t dy + Cjdz + ...+t l (3.1-3) 
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where 


x - x 0 ' dx, y = y 0 + dy, /. = z 0 + dz. ... 




^ b = Mi 
dx 1 dy ’ 



f t (Xq, yo, z 0 , . . . ) - L t . 


(3.1-4) 


Observation equation (3-1-3) can be written in the matrix form as 


v = Ax + 1. (3.1-5) 

That we have preferred to use weighted constraints to the station Chandler can 
be seen at the end of this section. These "a priori" weighted constraints on 
the station position generate observation equations of the form 


v* = Fx (3. 1-6) 

where F is a rectangular matrix, whose elements are either' zerojor one. Thus 
the complete observation equation system can be written as 


v = Ax + 1 
v x - Fx 


i. e. . 


where 


V = Bx + L, 



( 3 . 1 - 7 ) 


(3. l-7a) 


(3. 1-8) 


Observation equations for direction, Laplace azimuth, and distance are 
given in Appendix 11 fGrossmann, 1961, pp. 170, 177; Wolf, 1968, pp. 323- 
324], Due to angular and linear (distance) observations the observed data In 
a triangulation net is of a heterogeneous or dissimilar nature, 1 This hetero- 


lr Fhe term "heterogeneous or dissimilar" observations is used when the methods 
of their measurement are diverse; thus not only angles and distances, but also 
distances and heights are heterogeneous observations [Wolf, 1968, p. 561; 
[Schmid and J. Schmid, 1965 p. 10] uses the term "hybrid systems" for 1 'hetero- 
geneous systems". 
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geneous data have not only more than one dimension but also different 
"a priori" standard errors. To make this data homogeneous, i.e. , dimen- 
sionless and of unit weight, it is divided by the corresponding "a priori" 
standard error cr. For reasons of simplicity the mathematical model used is 
assumed to be uncorrelated. The resulting homogenized observation equation 
system can be written as 


v - Ax + 1 
v x = Fx 


(3.1-9) 


where 

v =- v/oi ; A = A/oi 1 = l/a 5 

v x = v x /q, ; Fj = F/oj (3.1-10) 

a t = standard error of L t ; qj = standard error of Xj . 


Equations (3. 1-9) and (3. 1-10) can be written in the matrix form as 


where 



(3. l-9a) 
(3. l-10a) 


respectively. 

Equation (3.1-9) is used directly for adjustment by conjugate gradients method 
in Sections 3. 23 and 3. 24 where these are used without ~ sign, although the A- 
matrix and 1 -vector used there are homogenized. A complete algorithm for 
obtaining solution vector and N 1 are given in Sections 3. 24 and 3. 26 respec- 
tively, which give v’Fv and Q xx or Q yy for a particular column. Using these 
quantities the "a posteriori" variance of unit weight (m^), standard errors (m x , 
m y )of unknowns, standard positional error (m p ) and the elements 0, A, B of the 
error ellipse are computed TGotthardt, 1968, pp. 121-125; Grossmann, 1961, 
pp. 163-168; Wolf, 1968, pp. 286-292]. Variance of unit weight m D 2 is given 
by 



m-n + c 


(3. 1-11) 
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Standard errors (m r ,'m y ) of unknowns are given by: 


m * - JQ X% ; m y = m 0 j Q yy , 


and the standard position error 


m 


/ 2 ~ 5 

= V m x + m y = m - 


d\ZQxx + Q 

The elements of the error ellipse are given by 


y y 


9 = — arc tan — 

2 Qxx ~ Qyy 


A = m 0 jQmx = semi-major axis of error ellipse 

B = m 0 J Q mln = semi-minor axis of error ellipse 


where 


(3. 1-12) 


(3. 1-13) 




Q, 


Bax, nln 


■- + ^1Y 


± 7 /(Q» 


Qyy) 2 + 4Q 2 


xy 


(3. 1 -14a) 

(3. l-14b) 
(3. 1-1 4c) 

(3. 1-1 4d) 


The standard errors of unknowns give the mean uncertainly of a station only 
in the direction of the coordinate axes, while the error ellipse gives this in 
any arbitrary direction. The standard positional error m p as well as the error 
ellipse possess an important characteristic that these are independent of the 
coordinate system rotation. 

The geodetic triangulation net is adjusted as an independent or free net, 
as it is not connected with other nets. For its unambiguoug determination, 
besides the observed data which includes directions, bases and astronomical 
observations, i. e. , longitude and azimuth, one fixed station is required to serve 
as origin [Gotthardt, 1968, p. 167; Grossmann, 1961, p. 175; Jordan/Eggert/ 
Kneissl, Bd. IV. 1, 1958, pp. 534-542]. Moses Lake station is kept as origin 
with its coordinates obtained from satellite triangulation results; these coor- 
dinates have been assumed to be the best known coordinates. 

Combining the free triangulation net with super-control net of zero order, 
i. e. , continental satellite net and/or super-transcontinental traverse means 
constraining the scale and/or orientation of the triangulation net. The effect of 
this combination is comparable with "tennis racket and string effect, " where 
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the rigid outer racket frame (super-control) constrains the loose strings 
(triangulation net). If the strings are already constrained, there would be no 
"visible" effect of the additional constrain from the rigid outer frame. This 
is also the purpose of this investigation, i. e. , to evaluate whether the existing 
geodetic triangulation is sufficiently "constrained" or needs to be constrained 
by additional super-control net. For the present investigation triangulation 
station Chandler, which is common to the three networks, provides constraint. 

Geodetic triangulation net can be combined with the super-control 
net in either of the two ways: 

(1) By using the actual data, i.e., by using the actual 
coordinates with their standard errors of Chandler as obtained 
from CSN and STT with the geodetic triangulation; or, 

(2) By adding a weight constraint to Chandler with its 
coordinates from the geodetic triangulation. 

For this investigation, the first way could not be used, as the super- 
control net coordinates of Chandler station are not compatible with those 
obtained from geodetic triangulation. As such, the second way has been 
preferred by using the actual preliminary accuracy estimates for Chandler, 
which are 1 part in 385,000 and 1 part in 3 million, as obtained from 
CSN ans STT, respectively. Further investigations are made by using 
hypothetical standard positional error accuracy estimates of Chandler 
station, which are 1:400,000; 1:500,000; 1:600, 000; 1:700, 000; ' 1:1 M; 

1:1,5 M. These accuracy estimates are within the actual preliminary 
accuracy estimates of super-control nets. Thus, using these various 
accuracies of super-control net, a feeling for the accuracy limit of super- 
control net, which would be necessary to improve the investigated geodetic 
triangulation, can be obtained. 
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3.2 Conjugate Gradient Method (Cg-Method). 

3, 20 Introduction. 

Although Cg-Method was developed by E. Stiefel and M. R. Hestenes 
independently from each other in 1952 f Stiefel, 1952, p. 23], this has been 
used only twice 1 ) for geodetic computations THilger and Remmer, 1967; 

Wolf, 1968, p. 1851. Its basic algorithm given in most publications THilger and 
Remmer, 1967; Schwarz, 1968 and 1970; Stiefel, 1952] is good for well- 
conditioned systems. In case of ill-conditioned systems (Section 3. 27) use of 
basic algorithm means too many iterations. However, use of certain 
formulas [Hestenes and Stiefel, 1952, p. 433] shows the suitability of 
Cg-Method even for very ill-conditioned systems. 

An attempt is made here to derive the Cg-Method and put its 
algorithm together because this is too scattered in mathematical liter- 
ature. The Method of Conjugate Gradients (Cg-Method) is a non- 
stationary relaxation method, which theoretically solves a system of 
simultaneous linear algebraic equations 

Nx + u = 0 (3.2-1) 

in n-iterative steps, where N is a symmetric and positive definite n x n 
coefficient matrix, x is a nx>l vector of unknowns and u is a n x 1 vector of 
constants. In geodetic work as the columns of m *n A-matrix (observation 
equation coefficient matrix) are independent (A-matrix is of full-rank), its 
quadratic matrix A T A(=N) is symmetric and positive definite. Then the 
system (3.2-1) - known in geodesy as the Normal Equations - has a unique 
solution. 

For the derivation of this method the matrix N of equation (3. 2-1) will 
be considered symmetric and positive definite, and then the derived algorithm 
will be modified for an arbitrary N-matrix and for directly using (he observa- 
tion equations without explicit formation of normal equations. 

1 ) so far known to this investigator. 
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3. 21 Basic Relaxation Method. 


Certain terms of the basic relaxation methods, which are used later, 
will be discussed first: 

To obtain the solution of (3. 2-1) by relaxation methods a trial vector x(°) 
instead of x in (3. 2-1) is used, which gives a residue vector r(°) 

r<°) = Nx<°> + u. (3.2-2) 

The aim of relaxation methods is to change the trial vector x(°) until the 
residue vector r^°) disappears. 

This residue vector r^°) is the gradient of the quadratic function F(x( 0 )) 
given by 

F(x<°>) = x<°) T Nx(°> + u T x(°). (3.2-3) 

- 

Differentiating partially the quadratic function F(x^ 0) ) w.r.t. the trial vector 
we obtain 

= Nx(°) + u, (3.2-4) 


aF(x(°)) 

ax(°) 


(3. 2-2) and (3. 2-4) give 

r (°) = grad F = Nx<°) + u. (3.2-5) 

Equation (3. 2-5) indicates that the solution of (3. 2-1) is synonymous with 
the problem of finding a minimum of the quadratic function F(x( 0 )) given by 
(3. 2-3) Tualston, 1965, p. 439; Schwarz, 1968, p. 45]. Let h be an arbitrar 
ily selected non-zero relaxation direction vector, which corrects the trial 
vector x(°) in the direction of h so as to achieve a minimum of the quadratic 
function F(x( 0 )), giving a new trial vector x', which is a linear function of 
the last trial vector x^°\ given by 

x' = x(°) + Ax(°> = x(°) + Xh (3. 2-6) 

where X is the relaxation distance factor or the correction factor for the 
unknown vector, which is so determined that the quadratic function F(x') 
be a minimum; F(x') is considered a quadratic function of the only variable 
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X with constant (fixed) x(°) and h. Thus 
F(x') = F(x(°) + Xh) 

Using (3.2-3) 

F(x') = ~ ((x(°) +Xh) T N(x(°) + >h)) i u T (x(°) *- Xh) 

2 

= — (x(°) t Nx^°) + >(x(°) T Nh) + ~ X 2 h T Nh + u T x(°) + Xu T h 
2 2 


Using (3.2-3) and (3.2-5) 

F(x') - F(x(°) + X r(°) T h + X 2 h T Nh. (3. 2-7) 

To find a minimum of F(x') in which the only variable is X, equation (3. 2-7) 
gives 

AE&-) = r (°) T h •+ Xh T Nh = 0 (3.2-8) 

Q A 


i. e.. 


r<°) T h 
h T Nh ' 


(3. 2-8a) 


Equation (3. 2-8a) means that the relaxation direction vector h cannot be 
selected orthogonal to the residue vector r(°), for in that case X Bln = 0 
at the trial vector x(°). The new residue vector r' is given by 


r' = Nx' + u 


(3. 2-9) 


It could easily be proved that the new residue vector r' must be orthogonal 
to the last relaxation direction vector h; in other words after each relaxation 
step the new residue vector r( J ) must be orthogonal to the last relaxation 
direction vector h^ J_1 ). 

Proof : 

From (3. 2-8) 

d F(x') • d F(x') dx' 

dX ~ dx' dX 
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(3. 2-9a) 

r( fl ) + XNh. (3. 2-9b) 

Hence 

r^h = (r< c ) + X Nh) T h = r<°) T h + X(Nh) T h 

r(°) T h 

= r(°) T h - TtT^T- (Nh) T h - 0 for X = \ nl „. 
h Nn 

i.e., 

r^h = 0. 

As proved above the new residue vector r' must always be orthogonal 
to the last relaxation direction vector h; this can be considered as an 
Orthogonality Condition . Equations (3.2-9a) and (3.2-9b) give the 
following relaxation distance factor X Bln : 

r n h = (r(°) + X Nh) T h = r(°) T h + Xh T Nh = 0 

_ r<°) T h 

mln h T Nh 

which is the same as given in equation (3. 2-8a). 

Going in the relaxation direction h with the correction factor X nln a 
minimum for the quadratic function F(x') is obtained; this can be proved by 
the second derivative of F(x') obtained from (8.2-8): 

d ^ = h T Nh > 0. (3.2-10) 

d X 

The second derivative in equation (3. 2-10) is always positive for every non- 
zero relaxation direction vector h, as N is positive definite TThomas, 1960, 
p. 130]. 


using (3. 2-4), (3. 2-5) and (3. 2-6) 

R. H.S. - r /T h - 0 

Alternative Proof : 

Using (3.2-5). (3.2-6) and (3.2-9) 

r - Nx / + u = N(x(°) + X h) + u = 


- 20 - 



Let AF be the decrease in the quadratic function F(x( 0 )). while going from 
x<°> to x in the relaxation direction h using X = X mln . Considering (3.2-7) 
and (3.2-8a) 

AF » F(xW) - F(x / ) = ~r\ r(°) T h + -^ X 2 h 7 Nhl = +7 > 0 * i 3 * 2 ' 11 ) 

' 2 2 h Nh 

for 

r*°) T h ft 0. 

(3.2-11) gives the largest decrease of F(x^°)) in the direction h. Thus by each 
relaxation step the current quadratic function F(x') decreases, which proves 
the convergence characteristic of the relaxation method. Geometrically the 
relaxation method can be interpreted as follows: 



Figure 3.2-1 


The quadratic function F(x( 0 )) in the case n = 2 can be represented by its 
level lines F(x( 0 )) = constant, which are concentric ellipses in rectangular 
(x u Xg) - coordinate system (Figure 3,2-1); the common center of these con- 
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centric ellipses coincides with the minimum point of F(x( 0 )), whose coor- 
dinates are (xj, x a ), which represents the solution of the system. Let 
x(°) be the initial point; according to (3.2-5) the corresponding residue 
vector r(°) will be orthogonal to the level line passing through x^°). A 
relaxation step means first to choose the relaxation direction vector h and 
then to proceed in the direction of h until x', where the quadratic function 
along the relaxation direction is minimum. This point x' is naturally on a 
level line where it is tangential to the relaxation direction vector. At this 
point x',the new residue vector r' is orthogonal to the level line i. e. to the 
relaxation direction h (3.2-9a). 

Based upon this basic relaxation method several relaxation methods have 
been developed, which differ in the selection of the relaxation direction vector 
h and the relaxation distance factor X for each step. 

3. 22 Method of Conjugate Gradients (Cg-Method) . 

Based upon the basic relaxation method, gradient methods were developed 
where the relaxation direction vector h is not chosen arbitrarily but is a 
function of the current or previous residue vectors. 

To find a minimum of the quadratic function F(x (°)) (Figure 3. 2-1) the quickest 
way obviously is to proceed orthogonally from x(°) in the direction h opposed 
to the grad F, i. e. , opposite to current residue vector r^ 0 ). This can be written 
in a mathematical form as: 

hW = -r( J_l) . (3.2-12) 

Equation (3.2-12) constitutes the principle of the Method of Steepest 
Descent. Although the decrease in the quadratic function in each relaxation 
step is locally maximum, the convergence in general is not good [Schwarz, 

1968, p. 68]. 

To improve the convergence, a method called the Method of Conjugate 
Gradients was developed by Stiefel and by Hestenes independent from each 
other, using not only the current residue vector but also the previous 
iteration results so that the system theoretically 1 has a solution in n-relaxa- 

x In practice due to round-off errors this theoretical convergence is not 
achieved (Section 3.40). 
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tion steps, i. e. it converges in n-finite steps 2 , n being the number of unknowns. 

Cg-Method is thus a modification of the Method of Steepest Descent, where 
the relaxation direction vector h( J ) is determined alone by the existing (current) 
residue vector r^ -1 ^ according to equation (3. 2-12), i. e. h( J ) = -r( J_1 ). In 
Cg-Method equation (8, 2-12) is valid only for the first relaxation step but the 
further relaxation directions are calculated by 

hO) = - r (J _1 ) + for j s 2 (3.2-13) 

( h( 3 ) = -r( J-1 ) for j = 1 i 

where is a Correction Factor for Relaxation Direction vector, "hich is 

proportional to the last relaxation direction vector h^ J_1 ). The factor e is 
so determined that the relaxation direction vectors h( J ) and h^ -1 ) are conjugate 3 
i. e. 

h(3) T Nh (j_1 ) - 0 (3.2-14) 

(3. 2-13) and (3. 2-14) give 

The relaxation distance factor Xj is determined from (3.2-8a) as: 

r( J_1 ) T h( J ) 

X s ~ ~ hO) T Nht J ) ’ (3 . 2— 1 (>) 

and the solution vector x (J) is determined from(3.2-G): 

x (0 = x 0“0 + h (j ). (3. 2-17) 

The residue vector r( J ) after j-th step is given by combining (3. 2-2), (3. 2-1 (!) 
and (3. 2-17): 

r< J ) = Nx( J ) + u = N(x^ rI )+X,h( J )) + u = (Nx( J-1 >+u) + XjNh( 3 ) - r( J ~ l > +X,Nh (j ) (3,2-18) 

2 This is why it is also called the "n-Step Method". 

3 Due to this relationship this method is called the Method of Conjugate 
Gradients. 
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The Cg-Method is thus defined by equations (3. 2-5) > (3. 2-12), (3. 2-13), (3. 2-15), 
(3, 2-16), (3, 2-17) and (3, 2-18) after choosing a trial vector x^ 0 ); generally 
x(°) = 0, 

Proceeding in the same way as was used to obtain equation (3.2-9a), it 
cari be easily proved that 

r (J)r h (J) = 0< (3,2-19) 

Similarly it can be proved that 

r (l)T h 0“0 = o. (3.2-20) 

Using (3. 2-13), (3. 2-14), (3. 2-19) and (3. 2-20) we obtain 

r (j)T r (j-i) = rO^-hM+e^hf 3-1 )) = -r^) T h^)+e J - i r( J ) T h^" 1 )= 0 (3.2-21) 

Equations (3. 2-14) and (3. 2-21) can be generalized in the following form by the 
method of induction: 

h^) T Nh^) = 0 fori ^ j (3, 2-l4a) 

r (0 T r 0) =o for i t j. (3.2-21a) 

From equations (3.2-14a) and (3. 2-21a), one can conclude that in the Cg- 
Method, the residue vectors r (3) constitute a mutually orthogonal system 
and the relaxation direction vector h (J) a mutually conjugate system. 

This gives an important characteristic of Cg-Method. As the residue 
vectors r (3) constitute a mutually orthogonal system in an n-dimensional 
vector space, this orthogonal system can contain a maximum of n non- 
zero vectors. Thus, at the latest, the (n + l)th residue vector r (n) must 
disappear, i.e. , r (n) = 0. This proves that the Cg-Method provides the 
solution in a maximum of n-steps. 

Numerical computations, however, show a deviation from the theoretical 
solution in n-steps, which may be due to (1) round-off errors and (2) an ill— 
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condition N-matrix; as these effect the mutual orthogonality of the residue 
vectors. Thus, the residue vector r (n) is not exactly zero. This deviation 
does not disturb the system. In this case where more than n- iterations 
are necessary to obtain the solution, the quadratic function decreases after 
every iteration. 

Equations (3. 2-15) and (3. 2-16) can be written in an easier form using 
(3. 2-13), (3. 2-18), (3.2-19) and (3.2-21): 

Numerator of (3. 2-16) = -r^ -1 ) 7 h( 3 ) = -r^ -1 ) T (-r( J-1 ) + e h^ -1 ) 

= + r (J-’ l ) T r( J - 1 )- €j _ 1 r( J “ l ) T hO" 1 ) - r( J_ 1 ) T r( J_1 ). 

Hence (3. 2-16) can be written as 

r( 3_l ) T r( J_1 ) 

‘ _ h( J ) T "Nh( J ) ' (3 ‘ 2_16a) 

From (3. 2-18) 

NhO" 1 ) = --i- (r^ 3-1 ) - r( 3-2 )) (3. 2-18a) 

Aj_i 

Numerator of (3.2-15) = r^ -1 ) T Nh^ _1 ) = — — — (r( J_1 ) - r^ -2 )) 

h-i 

= 7“ (rfl-OM 1-1 ) - r ( J_1 ) T r ( 3 " 2 )). 

«j-i 

r (3-l)T r (J-l) 

^i-i 

Combining (3. 2-15) and (3. 2-l6a) we obtain 

r (J-i)T r 0-i) 

£ 'J-i = 0-2). (3. 2-15a) 

Although the equations (3. 2-15) and (3. 2-16) are mathematically the same as 
equations (3. 2-15a) and (3. 2-16a) respectively, both sets of equations are 
equally good for well-conditioned systems. However, in case of ill-conditioned 
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systems equations (3.2-15) and (3.2-16) show better results [Hestenes and 
Stiefel, 1952, p. 433]. Consequently, these were the equations used in this 
investigation. Due to the simplicity of equations (3.2-15a) and (3.2-16a), 
they are the equations usually given in most of the publications. 

3.23 Modification of Cg- Method. 

1. If the A-matrix 1 is an arbitrary matrix and not assumed to be a positive 
definite symmetric matrix, Cg-Method can be used after multiplying the system 
by A T , i. e. 

(A T A)x + A T 1 = 0, (3. 2-22) 

for A T A is symmetric and positive definite, if A is non-singular matrix. The 
system of equations represented by (3.2-22) is equivalent to that given by (3.2-1). 

2. To use the Cg-Method for observation equations directly without forming 
normal equations (N-Matrix), i. e., without forming a symmetric positive 
definite system the following procedure is used. 

Let the observation equations be written as 

Ax +1 = v (3. 2-23) 

where A is a m x n matrix of observation equation coefficients, 1 is a 
m x 1 vector of absolute terms and v is a m x 1 vector of residuals. In 
case of observations with different weights, system (3.2-23) is homogenized 
by multiplying it with corresponding pj, where p* is the weight of the 
observations, thus the system becomes unitless and with unit weight. 

Normal equations to (3.2-23) can be written as 

A T Ax + a t 1 > 0 (3.2-24) 

which are equivalent to (3.2-22). 

Similar to equation (3.2-9) the residue vector r( J ) and the residual vector 
V' J ) after j-th iterations are 


^or reasons of consistency N and u of (3.2-1) has been replaced by A 
and 1, respectively. 
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r (0 = A T AxO) + A t 1 
and 

AxW + 1 = v( J ) 

Combining (3. 2-25) and (3. 2-26) we obtain 

{•0) — a t v^ 3 ^. 


(3. 2-25) 


(3. 2-26) 


(3. 2-27) 


In the algorithm of Cg-Method N can be replaced by A T A in equations 
(3.2-15), (3.2-16) and (3.2-18) and thus the expression 

h( J ) Nh( J ) = h( J ) A T Ah< 3 ) = (Ah( J )) T (AhO)). (3.2-28) 

Similar to equation (3.2-18) an expression for the residual vector v (J) after 
j-th iteration can be found: 

v( 3 ) = Ax( 3 )+ 1= A(x( J-1 ) + Xjh( 3 )) + 1 = (Ax( J_l )+ l) + XjAh( J ) = v( J_1 ) + Xj Ah( J ). 

(3. 2-29) 

Thus, the residual vector v (J 3 after j-th iteration can be computed in two 
ways according to equations (3.2-26) and (3.2-29), thus providing a 
computational check. However, to increase the stability of the relaxation 
process v (3) should be computed anew in each step according to its 
definition given by equation (3.2-26) and not recursively by equation 
(3.2-29). Thus, the residual vector rw is also to be computed 
anew according to its definition by equation (3.2-27) and not recursively 
by equation (3.2-18) [Lauchli, 1959, p. 259]. 

That the norm of v (3) decreases monotonously while solving the 
observation equations (3.2-23) by Cg-Method can be proved as follows 
from (3.2-29): 


y( J l ) = v( J ) - XjAh( J ) 
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(3. 2-30) 


vP _1 ) T v( J_1 ) = v0) T v0) -2Xjv( J )(AhM) + Xj a (Ah^^) T (Ah^)) 

From (3.2-27) and (3.2-19) 

v( J ) T (Ah( 3 )) = (A T v( J )) T h^) = r^ T h^ J V= 0 

(3. 2-30) can he now written as 

V 0)M J ) = - Xf(Ah( J )) T (Ah( J >) (3.2-31) 

Until x( J ~ l ) is not the solution x of (3.2-23), r( J ~ 1 ) / 0, hence € j — L j- 0 
and h(') i 0; thus Ah^ ^ 0 and Xj > 0, hence 

Xf(Ah< J )) T (Ah0)) > o. (3.2-32) 

Equation (3.2-31) together with (3.2-32) proves that the norm of residual 
vector v (3) decreases monotonously. 

3. 24 Algorithm of Cg-Method Using A-Matrix. 

Now a complete algorithm of Cg-Method for obtaining the solution vector (x) 
using directly the homogenized observation equations can be summarized in the 
following systematic way: 

Given : Homogenized observation equations Ax + l -- v. 


Select: Initial Trial Vector x<°) 

= 0 


Compute : 



(1) 

v<°> = Ax<°> + 1 




Relaxation steps j = 1, 2, . 

. . n 


(2) 

r( J_1 ) = a t v^ J ) 



(3) 

(Ar( 3-1 ) T (Ah( J "’ 1 h 

(AhO-^TjAh^- 1 )) 

( for 

j 2) 

(4) 

,(j) r-r(°) 

h -L-rP-O + e^hO- 1 ) 

(for 

(for 

i = i) 
j * 2) 

(5) 

, r< 3 ~ l >V J ) 

J (Ah( J )) r (Ah( j >) 
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(6) x( J ) = X ( J_1 > + \jh( J ) 

(7) v( J ) = Ax< J ) + i 

(7a) vO) ck = v( J_1 > + XjAhO) 
Tests; 

(8) Orthogonality Test: 

r( 3 > T h( J ) = 0 

r 0-i)T r (J) = o 



(9) = v( J ) 

' ' elm ok 


3.25 Termination of Iterations , 

Based upon the theory of Cg-Method and the geodetic requirements, 
iterations should be terminated as soon as any of the following conditions are 
fulfilled: 

(a) if’ the improvement: in the solution vector between two consecutive 
iterations is negligibily small, i. e., |x^)-x^ _1 ^| = 1.0- 10“ 4 seconds 
(i. e. 1. 0-10 4 second in <p or \ = 3.0mm), 


(b) 

(c) 

(d) 

(e) 


if = 0, 

if (Ah<*)) T (Ah(0) 0; 


if the given number of iterations is reached; 

if the round-off error (RFE) during.iterations exceeds a certain 

accuracy limit, which is given by the vector difference 


r - r 

1 true comp 


| 3 , where r = A t Axj + A T 1 = A t v( j ) and r Co 


A T v<’> . 

check 


RFE = |A T (vO)-v(O k )P 

The iterations should be terminated if r (3)T r (J) £3. RFE. 
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3. 26 N 1 - Inverse of Normal Matrix . 

To analyse a geodetic network variance-covariance matrix (S) is needed. 
This requires the computation of 1ST 1 , the inverse of the coefficient matrix of 
the normal equations, 

S = rho- Q where Q = N' 1 (3.2-33) 

By Cg-Method N 1 is computed column by column and for each column of N 1 
the computational steps needed are the same as those required to obtain the 
solution vector of (3. 2-1). For huge systems the computational time for 
obtaining one column of N x or to obtain the solution vector of (3. 2-1) is large 
and hence in such cases only a few necessary column vectors of N' 1 need 
to be computed. 

The algorithm for computing N' 1 is obtained in the following way: 

The inversion of N basically amounts to solving the system 

N*N*=E (3.2-34) 

or 

N- Q = E (3 . 2-34a) 

where E is a unit matrix. 

Let q k be the k-th column vector of Q (= 1ST 1 ) and the k-th column vector of E, 
then q k is the solution to the system 

N-q k -e k " 0. (3.2-35) 

q k , which is now the solution vector of (3. 2-35), can be obtained by equations 
(3.2-5), (3.2-12), (3.2-13), (3.2-15), (3. 2-16), (3. 2-17) and (3. 2-18) after 
choosing a trial vector q(°); generally = 0. 

Modifying this algorithm for using the A-matrix directly without explicitly 
forming the normal equations (N-matrix) N is replaced by A T A in the above 
mentioned equations. From (3.2-35): 

A T Aq k - e k = 0. 

Let qi°) be the initial trial vector, then 
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(3. 2-30) 


r(°) = A^q^ 3 ) - e k . 

The complete algorithm of Cg-Method for k-th column vector of N' 1 - inverse 
of normal matrix - using directly the homogenized observation equation 
coefficient matrix A can be given as follows: 

Given : Homogenized observation equation coefficient matrix A. 

Select : Initial trial vector q(?^ = 0 
Compute : 

(1) r(°) = -e k 

Relaxation Steps j = 1, 2, . . . n 


(2) 

Gj-i = 

(Ar^ -1 )) T (Ah( J-1 h 
(Ah( 3-1 )) T (Ah^ -1 )) 

(for j s 2) 

(3) 

h( J ) = 

-- r (°) 

(for j *= 1) 


__ r (J _1 ) + € J _ 1 h^ J_1 ) 

(for j s 2) 

(4) 

\ _ 

r^ -1 ) T h^) 



(AhW) T (Ah<0) 


(5) 

= 

q^ -1 * + Xjh^) 


(6) 

r (J) = r 0 -1 ) + X 3 A T (Ah^)) 



Test : 

(7) Orthogonality Test: 

Same as (8) of Section 3. 24. 

Termination of Iterations: Same as in Section 3. 25. 

3. 27 III- Conditioning . 

If small errors in the coefficients of equation (3. 2-1), i. e., equation (3. 2-24) 
or in the solving process have little or no effect on the solution the system is 
called well-conditioned; if the effect is large it is called an ill-conditioned 
system. Ill-conditioned systems have a very poor rate of convergence. A 
system can be evaluated if some information about its condition (condition 
number) is known. Condition numbers can be computed by using equations 
given in [Fox, 1965, p. 142; Schwarz, 1968, pp. 22-23; Zurmuhl, 1964, 
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pp. 212-214], For the present investigation condition number (K) is the ratio 
of the largest eigenvalue (Xna X ) of N to its smallest eigenvalue (X slB ), i. e., 

K = (3.2-37) 

This condition number K should be near unity for a very well-conditioned 
system from the point of view of solving linear equations [Fox, 1965, p. 199]. 

As mentioned in Section 3.1 eigenvalues need not to be calculated for 
determining the convergence as Cg-Method is a finite iterative method. This 
characteristic (i.e.,no computation of eigenvalues) is of special importance 
in the case of modified Cg-Method as the original A-matrix is directly used, 
where due to the lack of a square matrix (N-matrix) computation of eigen- 
values is not possible. In spite of this a condition measure can be derived 
bv using the Cg-Method algorithm in the following way: 

In THestenes and Stiefel, 1952, pp. 416-420] it is shown that 

T~ < X i <T- (3.2-38) 

''■DBX ^Bln 




(3. 2-39) 


T Nh( 3+1 ) _ _Xj_ h^~ l ) T Nh ( ’> 

h( J ) T Nh(^ “ Xj-l’ h( J "‘) T Nh< 3-1 ) 



The last two equations show that if, at the beginning of any iteration, 
r (j - 2 ) t r (j - 1 ) o an( j h(.i-i) T NhO) ^ 0 then the computed values of r( J-1 ) T r( J ) 

and h( J > T Nh( J+1 ) will also deviate from zero. This deviation will depend upon 
the magnitude of ALPHA = X^/Xj-f, the larger this ratio is the greater will 
be the disturbance of the orthogonality relations, and more rapidly the round- 
ing-off errors accumulate. According to equation (3. 2-38) the condition 
number K (= Xaax/^m\ n ) is an upper bound of the critical ratio X 3 /Xj_ 1( which 
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determines the stability of the process. When K is near one, that is, when 
N is near a multiple of the identity, the Cg-Method is relatively stable 
f Beckman, 1960, p. 70l. 


3. 28 Round-Off Errors . 

Round-off errors could be mainly due to two reasons: 

(1) For adjustment, the approximate values of the station coordinates are 
to be used in decimal degrees, where these decimal values are rounded-off 

to a certain significant digit; 

(2) Due to the limited number of digits per storage location in computers, 
round-off error occurs during arithmetic operations. 

A detailed study on round-off errors is given in TMuller-Merbach, 1970], 
wherein the investigations show that only double precision storage combined 
with double precision computations will keep the round-off errors to a 
minimum. To minimize the computer-error while adding small numbers to 
large numbers, the summation of scalar products r^ 3 ) T r^), v T Pv and 
(Ah( 3 V(Ah( J )), which are used for iterations, is done in blocks. Investigations 
on a system of 66 equations and 39 unknowns show that r (3)T r (3) and (Ah (J) ) T 
(Ah^ J) ) are very sensitive to the block size while v T Pv is quite insensitive to 
it (Table 3.2-1). It has been found that the best convergence, that is, 
solution vector for r( J ) T r( J )=0 after least iterations, is obtained when the ratio of 
blocksize (KMM) for (Ah( 3 )) T (Ah( J )> to that (KM) for r( 3 ) T r( 3 ) is afunction of the ratio of 
the number of equations to that of unknowns. This has been programmed in 
Subroutine PATSUM. The only parameter to be determined in this subroutine 
for individual problem is the basic blocksize (KM) for r( 3 ) T r( J ). 

A characteristic measure for the round-off error (RFE), which is required 
for the termination of the iterations, is the norm TGinsburg, 1963, pp. 197-198] 

RFE = |r tru , - roonj 3 (3. 2-41) 
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where r tri „ is the true residual vector by definition according to equation (3. 2-9) 
or (3. 2-27) and r coBl , the recursively computed residual vector according to 
equation (3. 2-18). Hence 

RFE = |A T (v^)-v( c J 2 ck )| a (3.2-41a) 

which provides a condition for termination of iterations, given by 

r( J ) T r< 3 ) < 3. RFE. (3.2-42) 


Table 3.2-1 

Effect of Block Size on r (3)T r (j \ (Ah (J >) T (Ah (J) ) and 
V T PV at 80th Iteration 


C2 

g 

Block 

Size 

r 

i 

i 


•i“-d 

Si 

CO 

a 

X 

W 

KM for 
RTR 

KMM for 
(Ah (J V(Ah (j) ) 
or V T PV 

r O) T r (j) 

(Ah a) ) T (Ah (j) ) 

V T PV 

a 

10 

10 

0.4127 D03 

0.2020 D10 

0.2529 D01 

2 

10 

15 

0.1594 D02 

0.7879 D08 

0.2529 D01 

3 

10 

20 

0. 1087 D03 

0. 1286 D10 

0.2529 D01 

4 

10 

NE 

KM ^ 

0. 1247 D03 

0.6080 D09 

0.2529 D01 


NE = Number of Equations 
NU = Number of Unknowns 


3 . 29 Residue Vector (r) . 

As mentioned earlier each approximation to the solution vector is 
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closer to the true solution vector x than the preceeding one, the error vector 
(x-x( 3 )) also decreases at every step. However, the residue vector r may 
not decrease; normally the "residue square" |r| s oscillates and may even 
increase. This oscillating behaviour of the residue square lr| s could be 
due to round-off errors. Investigations have shown that the solution vector 
xO) is no longer improved substantially as soon as the norm |r| 2 comes down 
to the range of the round-off errors. This occurs when the iterations are 
continued beyond the number of unknowns, i. e., j >n, when the residue vector 
r^) will sooner or later begin to decrease sharply such that ultimately x( J ) 
will be as good an approximation to x as it may be expected from the condition 
of the system [Beckman, 1960, p. 69; Ginsburg, 1963, pp. 192.199; 

Hestenes and Stiefel, 1952, p. 411]. The following systems were investigated: 

System I: 39 unknowns, 66 equations (Figure 3.2-2) 

System II: 84 unknowns, 138 equations (Figure 3.2-3) 

System III: 573 unknowns, 965 equations (Figure 3.2-4) 

In Figures 3.2-2, 8. 2-3 and 3.2-4 asterisk marks (*) show that the current 

solution vector has no substantial improvement 'over the last solution vector, which 
happens after such n-iterations where the round-off errors stop increasing. 
Thus, any of these solutions are as good an approximation to the final 
solution vector as could be expected from the condition of the system. 

3.3 Programming 

All programs were written and tested by the investigator, except the 
SUBROUTINE DUMMY, by F. Fajemirokun. 

3.30 A-Matrix, Index Matrix and Storage Requirement 

Structure of A-Matrix and F-Matrix. 

The homogenized A-Matrix, in general, contains m rows and n column s 
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where m = number of observations and n = number of unknowns, which is 
equal to 3 times number of stations, as for each station there are three 
unknowns, namely orientation correction (dz) and the two coordinate corrections 
((Ur and dX). Each row of A-Matrix contains a maximum of five non-zero 
elements, which are so arranged that the first term is the coefficient of orienta- 
tion correction, the other four terms are the coefficients of station coordinates 
corrections. Each row of F-Matrix contains a maximum of one non- zero element. 
Thus for a large geodetic triangulation net there will be much more zero elements 
in each row than the non-zero elements, i. e. , A-Matrix is a sparse matrix. 

The observation equations for the same station are put together one after another. 
Index Matrix and Storage Requirement . 

The original A-matrix is a sparse matrix and contains a maximum of 
five non-zero elements in its each row. Thus in a large triangulation system 
of 1400 observations and 820 unknowns the non-zero elements are only 0. 6% 
of the original A-matrix. To save the computer storage space for storing 
the original A-matrix, Index Matrix (II) has been used so that the "reduced" 
A-matrix has only five columns. The elements of the Index Matrix are positive 
integers. A very simple algebraic expression was deduced to obtain the 
elements of Index Matrix from the A-matrix; the lowest element of Il-matrix 
is equal to the lowest station number and its largest element is equal to 3 times 
the largest station number. A program gives?!! -matrix from A-matrix. 

Thus Index Matrix, also called' Reference or .Guide- Matrix^ reduces 
an original (m x n) A-matrix to a "reduced" (m x 5) A-matrix. Both 
the "reduced" A-matrix and the Il-matrix together need much less storage 
than the original A-matrix. For example for an original (1400 x 820) A-matrix 
114800 storage locations and for N-matrix 336610 storage locations stored in com 
pressed form are needed, whereas the "reduced" A-matrix together with Il- 
matrix needs 14000 storage locations, which is 1.2% of the original A- 
matrix and 4. 2% of the N-matrix. This would mean that the same computer, 
which cannot handle an original A-matrix or N-matrix, can. easily handle 
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the same system using Index Matrix, as well as can also handle much 
larger systems. In future discussion, the term A-matrix will be used 
for "reduced" A-matrix. 

3.31 Minimizing Round-Off Error. 

The ideas of Section 3. 28 were used to minimize the round-off errors. 

Double precision storage was used together with double precision computation; 
for integers single precision was used. Summation of scalar products 
r(') r r(''), (Ah^V(Ah^) andv T Pv was done by using the subroutine PATSUM, 
which minimizes the effect of adding small numbers to large numbers. 

3. 32 Solution Vector Program and N 1 - Program . 

The algorithm of Sections 3. 24 and 3. 25 was programmed as a SUBROUTINE 
SOLN. The algorithm of Section 3. 26 was programmed as a SUBROUTINE QSOLN. 
These programs were initially tested on five data sets. The results obtained by 
these two subroutines were found to be the same (within computational accuracy) 
as those obtained by using Gauss-Algorithm. The five data sets used were: 

(3, 2), (10,4), (16, 12), (36, 24) and (66, 39) - the first number = number of 
observation equations and the second number = number of unknowns. 

The solution vector or N^-column will be printed after every N-iterations 
(N^Interval for writing computed output) besides after the fulfillment of the 
usual conditions. 

3. 33 Universality of Programs . 

Both subroutines, namely SUBROUTINE SOLN and SUBROUTINE QSOLN, 
can be used for any feasible size of data, which can be accomodated on the 
available computer, after changing KM, which is the PATSUM Basic 
Block Size for RTR. 

The main program used together with these subroutines has dimension 
statements and a data card for Number of Unknowns (NU), Number of 
Equations (NE) and Number of Columns of Index Matrix (NI), which can 
be changed if there is need for it. 

The program is universal in the sense that it can be used for varying 
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data without much change and that "mesh-point numbering technique" is 
not required. Therefore, stations can be added or taken out from the 
triangulation system without worrying about the band-width and size of 
blocks. These programs have been tested on systems from as small 
as 2 unknowns, 2 equations up to as large as 804 unknowns, 1397 
equations . 

3.4 Number of Iterations and Computer Time 
3.40 Number of Iterations. 

As mentioned earlier in Sections 3.1 and 3.22, the Cg- Method theoretically 

gives the solution vector at n-iterative steps (n = number of unknowns) 

[LSuschli, 1959, pp. 257; Schwarz, 1968, pp. 74 and 1970 pp. 133; Wolf, 1968, 

pp. 184]. Investigations, however, show that the solution vector is not 

achieved in n-iterations due to round-off errors, ill-conditioning of the system, 

disturbances of the orthogonality and of the conjugacy relations. In [Beckman, 

1960, pp. 69; Fox, 1965, pp. 210; Hestenes and Stiefel, 1952, pp. 411] it is 

mentioned that frequently (n + l)th solution vector is significantly better than 

the nth one. According to [Ginsburg, 1963, pp. 192] up to 3n iterations may 

be needed in case of bad conditioned system; while [Hilger and Remmer, 1967, 

3 

pp. 13-14] mentions that n to — n iterations are needed in case of large systems 
for 4 decimal accurate solution vector and that number of required iterations 
strongly depends upon the condition of the matrix A T A. 

This investigation, using the actual data set (Section 2), shows that the 
number of iterations required to obtain the solution vector by Cg-Method 
using directly the A-matrix without explicitly forming the N-matrix depends 
upon two factors: (1) condition of the system, and (2) accuracy of the solution 
vector required. The dependence of the number of iterations on the condition 
of the system has already been discussed in Section- 3.27. Using 
the geodetic triangulation data (573 unknowns, 963 equations), the program 
went up to 5778 iterations without giving any 7 decimal accurate solution 
vector, while 4 decimal accurate solution vector was obtained after 1161 


- 43 - 



iterations, i.e. , 2.1 times number of unknowns (Table 3 .40-1). 

Each column vector q k of N' 1 is generally computed in less than 
1;2 n-iterations (Table 3.40-1). 

Table 3.40-1 


Experiment 

Number* 

Number of 

Solution Vector 

Covariance Vector 
for Column 8 

Unknowns 

Equations 
■ ■ ■ ■ 

Iterations 

Time** 
m sec 

Iterations 

Time** 
m sec 

1 

573 

— — — — - 
963 

1161 

9 

37.13 

640 

3 

45.96 

2 

573 

965 

1177 

9 

23.27 

657 

3 

31.91 

3 

573 

965 

1175 

5 

45.97+ 

659 

2 

12.59 + 

4 

573 

965 

1176 

9 

22.32 

682 

3 

45.64 

5 

573 

965 

1164 

5 

53.44+ 

674 

2 

1.77+- 

6 

573 

965 

1162 

5 

41.16+1 

675 

2 

0.00+ 

7 

573 

965 

1166 

9 


631 

3 

20.03 

8 

573 

965 

1159 

9 

24.29 

648 

3 

19.29 

9 

573 

965 

1169 

9 

29.41 

608 

3 

11.51 


*Refer to Table 3.5-1 

**Time is the Execution time on H-Compiler, Option = 2 (IBM 360/75) except 
those marked with a plus (+) sign, which is the Execution time on H-Compiler, 
Option = 0 (IBM 370/165). 

3.41 Computer Time. 

As for this investigation, IBM 360/75 and 370/165 are used. Factors 
influencing the computer time are valid only for these types of computers. 

For large systems, H-Compiler with Option = 2 (IBM 360/75) is found 
to be approximately 10 times faster in execution than the G- Compiler, 
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although compilation takes a longer time in H- Compiler than in G- Compiler. 
For (819 x 1395) system, G-Compiler took 9.58 minutes for 98 iterations, 
while the 11-Compiler took 7.00 minutes for 791 iterations. The time 
given in Table 3.40-1 is thus for 11-Compiler. It is noticeable that time 
required for ohe column covariance vector is much less than the time 
required for solution vector. 

3.5 Results. 

The results of the investigation are given in Tables 3.5-1 and 3.5-2, 
wherein the improvement of the particular geodetic triangulation by super- 
control net is visible only when its accuracy is at least 1 part in 500,000. 
The positional improvements of Wyola (95), which is in the middle of the 
triangulation chain, using various station constraints for Chandler (3) 
are shown in Figure 3. 5-1. These positional improvements are relative to 
free net adjustment results. As the preliminary accuracy of continental 
satellite net is lower than 1 part in 500,000, this cannot be useful as a 
"constraint" to the geodetic triangulation net. On the other hand, the high 
accuracy of super-transcontinental traverse, which is one part in 3 million, 
makes it most suitable as a "constraint" to the geodetic tri angulation net. 

Worth mentioning is that the longitude terms, which are Q yy and m y s 
in Table 3. 5-1, remain practically uneffected during the entire investigation. 
This could be explained by the fact that station Wyola is very close to Laplace 
stations, which control the azimuth error accumulation, thus effecting the 
longitude error TBomford, 1965, pp. 90, 519]. Hence, due to closeness of 
Laplace stations, the longitude terms remain practically uneffected. 
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Table 3. 5-2 


c 



WYOLA (95) 


g © 

© s 

Q, 3 

Accuracy 




Positional Improvement 
Relative to Experiment 1 

X 55 
W 

1 in 

m x 

m y 

m p 

Meters 

% 

B 

Free Net 

1.83 

0.37 

1.9 



2 

300,000 

1.93 

0.37 

2.0 

-0.1 

- 5 

3 

400,000 

1.81 

0.37 

1.8 

0.1 

5 

4 

500,000 

1.51 

0.37 

1.5 

0.4 

21 

5 

600,000 

1.51 

0.37 

1.5 

0.4 

21 

6 

700,000 

1.51 

0.37 

1.5 

0.4 

21 

7 

1,000,000 

1.43 

0.37 

1.5 

0.4 

21 

8 

1,500,000 

1.33 

0.37 

1.4 

0.5 

26 

9 

3,000,000 

1.08 

0.37 

1. 1 

0.8 

■ 

42 


Standard Errors of Unknowns (m x ,m y ) and Standard Positional 
Error (m p ) are given in meters. 
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4 . SUM M A R Y A ND C O NC LU SIO NS 


The super-control net, i. e. , continental satellite net or super-trans- 
continental traverse, can provide a useful constraint to the investigated 
geodetic triangulation net, and thus can improve it only when the accuracy 
of super-control net is at least 1 part in 500,000; in this case, this 
corresponds to ±3. 7m standard position error for the station Chandler. 

The preliminary accuracy of super-transcontinental traverse is already 
better than this limiting accuracy of 1 part in 500,000. The preliminary 
accuracy of continental satellite net is, however, lower than the limiting 
accuracy of 1:500,000; the preliminary standard position error for Chandler 
as obtained from continental satellite net corresponds to ±4. 8m, i. e. , 
1:385,000. The futu re will show whether the limiting accuracy could be 
achieved by continental satellite net, especially because numerous spatial 
triangulations of CSN have produced accuracies within the range of 1 part 
of 400,000 and 1 part in 700,000 f Schmid, 1965, p. 22]. 

Schmid I"] 970, pp. 23-24] indicates that continental satellite net will fall 
short of an optimum solution with respect to both it's coverage and its 
accuracy. The three-dimensional positions of CSN-stations will probably be 
determined to no better than ±4 meters in all components, which does not 
seem to be good enough at least for this particular investigation. 

It might be useful to have a "block constrain" instead of "chain constrain", 
that is, to use four well separated satellite stations, namely 003, 102, 112 and 
1.34 (Figure 2. 0-3) over a very large area, thus constraining the triangulation 
of the western-half of the United States instead of one triangulation chain 
("chain constrain") between stations 003 and 102. 

Super-transcontinental traverse can provide a better constraint, if more 
then two of its stations are common to the stations of geodetic triangulation net. 
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Also, a "block constrain", as explained above, might be more useful instead 
of a "chain constrain". 

The development tendencies of instrumentation indicate that the future 
super-control nets will use VLBI (Very Long Baseline Interferometry) and 
Laser ranging systems. 



APPENDIX I 


Helmert-Rainsford-Sodano Inverse Problem Solution 

As the distance between stations varies' from 273m to 190km Sodano's 
concept which gives formulas applicable to very short as well as very long 
lines to solve the Inverse Problem has been used. The formula renders 
accuracies on the order of Of 00001 in azimuths and a millimeter in distance 
for any length of lines fHelmert, 1880,Vol. 1; Jordan/Eggert/Kneissl, 1959, IV. 2; 
Rainsford, 1955; Rapp, 1969; Sodano, 1958, 1963 ]. The iteration process 
stops when the value of (X- L) (i.e. , longitude on the reduced sphere - longi- 
tide on the ellipsoid) does not differ by its preceeding value by Of 00001. 

GIVEN: 

B,, L t Geodetic latitude and longitude of a point P x 
Bg, Lg = Geodetic latitude and longitude of another point P 3 

(longitude positive eastward, latitude positive northward). 

REQUIRED: 

Geodetic distance S 12 , direct azimuth A IS and reverse azimuth A al 
(Azimuth reckoned clockwise from north). 

NOTATIONS USED: 

L = Lg - L, = Difference of longitude on the Ellipsoid 
X = Difference of longitude on the reduced sphere, for which a 
progressively better value is calculated after each iteration 
(for first approximation X = L) 
a = Azimuth of the specific geodesic at the equator. 

CALCULATIONS: 

1 . Reduced latitude for each point: 

= arc tan r (l - f) tan Bjl. 
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2. Spherical arc 0 between two points 

sin a = r(sinXcos jS 2 ) s + (sin 8 2 cos|8i - sin#! cos # 2 cosX) 2 ] s 
For first approximation let X = L = Lg - Lq 


3. 


sin X, cos fj, cosS 2 

sin » = L - : — cu 

sin cx 


cos 2 or B = cos cr - 3 


2 sin 8 1 sin/Sp 

COS 3 (7 


From this relation further relations for cos 4a m , cos 6cr m and cos8cr B 
can be deduced: 

cos4cr B ,= 2 cos 2 2ct b -1; cos6ct b = 4 cos 3 2 ct b - 3 cos 2 ct b ; 
cos8ct. b = 2 cos 2 4ct b - 1. 

5. (X- L) = f sina(Ao<j + A 2 sincr cos 2ct b + A 4 sin2cx cos4ct b + 

Ae sin 3cr cos 6cr B + • • • ) 

where 

A o = l-'J ffl+f + f 2 ) cos 2 ry+ f (1 + | f)cos 4 a-^ f 3 cos 6 O' 

1 IQ 7^ 

A 2 - — f (1 + f + f 2 ) cos 2 o' - — f 2 (1 + — f) cos 4 a + — f cos 6 a 

4 4 4 256 


A ‘ * H fS < 1+ T 00084 a "xk 


f 3 cos 6 tt 


^6 


768 




COS 'O'. 


6. Iterate till (X- L) does not differ with its last computed value by 1. 10 12 
radians (as double-precision computation is done). 

7. X = L + (X - L) 

sin X cos 

8. tan A i2 = — : — — 1 . — — 

sin 8 2 cos Pi - cos X sm cos S 2 


tan A 2 i = — : 


sin X cos fjj 


sin S 2 cos fii cos X - sin j3 l cos $ 2 
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The following table is used in the assignment of the quadrants for azimuths 


Sign of L 
L = L 3 - L, 

Sign of 
tan A [2 

Quadrant of 

Aia 

Sign of 
tan A ai 

Quadrant of 

A-21 

+ 

-t 

I 

+ 

III 

+ 

- 

II 

- 

tv 

- 

+ 

m 

+ 

I 

- 

- 

IV 

- 

II 


9. S is = b(B 0 cr + Bj sincr cos 2cr m + B 4 sin 2 cr cos 4 cr n + sin3cr cos 6a n 

+ Bg sin 4cr cos 8a m + • • * ) 


where 


Bo = 

j. 

1 4 

U 2 _ u 4 + 

64 

— u 6 - 
256 

16384 

B a = 

1_ 

,.2 X ..4 

_15_ „6 + 

25 

4 

16 

512 U + 

2048 ’ 

b 4 = 


J- + 

3 6 

35 


128 

512 

8192 


.8 


Be = 


1536 6144 


u 


Be 


65536 


where 

u 2 = K 2 = e' 2 cos 2 oi 
b = semi-minor axis of the ellipsoid. 


This completes the Inverse Problem solution. 
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APPENDIX II 

Observation Equations (for Ellipsoidal Geodetic Coordinate System) . 
NOTATIONS USED: 

Bj , Lj = adjusted coordinates of station j (where j = i or k) 

B°, L° = approximate coordinates of station j 

Mj , Nj = meridional - , normal radius of curvature at station j 

Z, =. adjusted orientation unknown at station j 

L lk = observed direction from station i to station k 

\*j = observed astronomical azimuth and longitude at station i 
s lk = measured distance of the geodesic ik 
A° k = computed ellipsoidal azimuth of the geodesic ik 

s° k = computed distance of the geodesic ik from approximate coordinates of 
stations i and k 

Z° = computed value of orientation unknown Z i by "mean-orientation-method" 
dBj, dLj = corrections to approximate coordinates of station j 
dZj = correction to computed orientation unknown Z° at station j 

v ) = number of observed directions at station j 
«t lk = absolute term of the observation equation 

v lk = residual to the observed quantity. 
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Direction Observation Equation fGrossmahn, 1961, p.177; Wolf, 1968, pp:275, 323] 
v ik(pj r ) = ~ dZ i + “I 1 - sinA° lk - dB 4 - (^ . cos • cos B° - sinB^-dLj 


’lie 


5 lk 


+ ^ • sin A^! • dB^ - ■ cos A^ • cos B^ • dl^ + l. 

s ik s lk 


lk 


(II. 1) 


where 


^ik ~ A° k - (L lk + Z?) 


yO _ _L r A 0 _ r -l 
/j l ' rt lk ^lk I 

V\ 


(II. 2) 


Laplace Azimuth Observation Equation TGrossmann, 1961, p. 177; Wolf, 1968, p.323]. 
v iw LaD ) = sinA°.dB! - % ■ cosAf k -cos B?-dLj + ^ • sin A& • dBfc 

\ua p 7 s lk »ik S lk 


. K. 


Jo • cos A^t • cos B^ ■ dL k + l i]t 

s lk 


(II. 3) 


where 


= (A£ - otf) - (L° - X?) sin B? (II. 4) 

Distance Observation Equation TGrossmann, 1961, p. 170; Wolf, 1968, p.324l. 


r ik( D18 ) -^•cosA° k ‘dB 1 -^ L *sinA° k *cosB°-dL 1 - cos A?,- d^ 


-^-•sinA^-cosE^-dl* + (II. 5) 


where 


Ljk — s 


ik ^ik 


(II. 6) 


The corrections dB^ ,dL, and dZ , are obtained by the method of least squares; thus 
the adjusted station coordinates (Bj , Lj) and orientation unknown Zj are 
given by 
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(II. 7) 


Z J = Z?-.-dZ J 
Bj - Bj° + d R, 

L, = L° t dl, 

DIMENSIONS: 

All angular quantities are in seconds of arc and linear quantities in meters. 
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APPENDIX III 

Programs for Solution Vector and N' 1 by Cg- Method 
Due to comment cards, the programs are self-explanatory. 
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c 


c 


L(<?qn),M(O 90 ) , 


C 


HA IN* PROGPAP. OP r.G-MFTHnn F 6 R C OVA R ; I A N C F MATR IX ' 

I” Pi. TC IT R FAL * 9 ( A-t! f n-7 ) ' ' MyM 

P 1 MFUS I nr.; a (qc: 0 , M,I).( c 90, 5 1, X ( A.0.0 ) , c ( AGO )., 

1 LK ( 9*0 ) , MR ( Ciqn ) , I L ( <?°G ) , KL ( <?'O'0 ) , R ( 600 ) , H ( AGG’I , AH ( 990 ) , P, T ( 1 , AGO )', 

? R I. S T. ( A 00 ) * XL c T ( ago ) , ar.( c.Qf> ) - 

VAR I AOLFS nr FA'S rn.V ST ATEMFNT fM ;THF ’ A! At M*"Vo O.nCR A M * 

V /' I! I A f I. F : DATA OAR 0 FOR Mt.l : OF >1 

Mil mijmp.fr OF nr'KMOi.if.i<; ■ 

I l ; - M'IMPRI? no np$FR VAT I OM FOMAT TOMS 

M I - Ail IMP FR OF COLUMNS IM'IMOEX MATRIX ( tl V " “i 

R FAD ( F i 5000 ) ,viT,me,MJ M V ■ 

F ooo format <?ifv" '* rr- - - •- - v ~~ — ? 

' CALL OROLMI A , ) 1 , X ,F , M L , M , LK , MK . I L » KL , R , H , AH , R T , R L ST , XL ST, 

1AR ,'- H, MF,ni ) -- 

STOP _ 

RHP POUT IMF ■ O.SOI.M ( A ,T I /X , c * ,L >T-'» LKVA'''< , T L , KL , R » H V.AH , RT ,PLST,XIST',' 

iar,hu,ivf,mj j .'• - ■' '• ; -'ll.:.. - ;, ; 


0 

c 

c 


0 

c 


COVARIANCE MATRIX RV CO— METHOD 
IMPLICIT PFAt.*0( A-H,0-7 ) ’ V 

■ HIMFHRinr' Af* ,r -,MI V, n (HFVMiriX'fMII )', F(,liU), - --- LCMF.IVKfMF) , 

] L K ( !-' F ) t M K I ME )» !L (N'F ) ,KL:U‘F)*R (MW) , H I Ml J > ♦ AH ( NF|,RT ( 1» NU ) » R I. ST ( MU ) 

?, yLRTI Mil) , AFMMF ) -r--- - - - -y ' : “v 

VARIABLES IM S HP R OHT IMF- . fi,M KM ’ ;••• 

H I FF = l.o n-07 - 

M-v - ?*i s ip . .. . 

mi- 6 .= 6*Mf ” ;~ !T "p •* 

fmjn = l.o o-o° : . .c. ,, ,V 

M n IMT.FRVA'L FOP HP ITF-nilT PUT" (STATEMENT' MRill A) 

N - AOO V -;it 

Kii = CODE FOR A D JUSt F 0 RT AT I OM" CnORD . PllHOHEDyOUT PUT 

kp = <‘5 . : . 'Mo. 


Kf - PAT SUM BASIC P L OCK S I ZF < FOR^RTR ) 

FKTP. = DFLOA T ( MF ) /OFLOAT (AH* ).><> 

' K'-’li =IOIMT ( FKTP*nFLOATCKM') V X-syt" ' " V" 

1 1- 1 PI IT : A -MAT R 1 X At-'n II -MATRt.Xl/l;. ; 

RFi.'TMO l "" . . \ • 

IT 11 J = 1»mf /v .• 

p . f a n ( i > - fTJ' ) -» Tt 1^:?^ y k = V t*C C .U ^ 

If K( J ) , ( I 1 ( J,K ) ,K=1 ,5 ) ,K’L iil ' ) • ■ 

1 1 C nM T I r ■ i IF " ^rr^-T-T “ -• -rv-M-T-'M-; •- 

or* in j - i ,nf 

if ( i.m ) .mf.lki j v.np.Hf jrrnjr. il ( ji .mfVxl un eo to ?soo 

]0 COAT JMJF ... • y. . .. '..M ; 

HR IT F ( fi , AOOO I ( L C J J , r* f J ) , ( ACJ, KJ,K=1, F) , IL( J)>L k ' < >’ J . ' 

If- K ( J ) , ( ! 1 ( .J. K I tK“l » r '1 »KL ( J )VJ=T*MF) ■ 

UP I TF ( 6 , A,nni ) M H , OF - - 

CALL CSC A. L F ( l.o' r-np , a , A'F t m i ) 

1 p, c a. f) ( fj.» 5 ?on.) ''Ffin ; Trrrry^T^rr^r-’-.^ - — . - — 

lFIRFi--11.Er.7n9c) r-n to 9000 . ■ t u r > 

wpitf(a,'oi?) wk :'-o : 

fin 12 J = I , Mir: * . ---r- — r-*y*- - - 

1? F(j) = 0.0 no . . - ..... . ■';/ •; : /KJ/ ■ 

fio 16 J. = 1 , 0(1 . '< •?': .. ■ 58 ‘ ; ;. 

- R( j ) = o.n no - ... 



ri.ft ( j) = n.c no 
16 x ( j ) '=' ;o:o"do 

CALL MSCALF l 1. .0 0-04 t E t K'U* 1 ) •. 

00 15 • J i 1 « I'M I - - -- 

IF R C J ) = 

flolst = o.o no -r 
KPIJNT - ' ' • 

KKOUMt = 1 ; " 

RTR .= o.o no 

RTP.LST = O . o DO - 


OP 4 6 K a 1 , Ni l 
"4" 6" XL ST CK )' = 0'. 0;"00 

CALL PA TS t ( R , Ml ) , KM , RTR ) 
WRITH 6,50? ) RTR ,KOIIi\ ! t~"~” 
]F( RTRLST .F O.O.o) op TO PI 
4 5 FLA ST = ARTAH/AHTAH T" 


PI RTRLST a RTR 

IF(KOUNT-n 85,6 5,8.8 : — ‘ 

F 5 DP 86 I = ], Ml l 

86 H(I) = -R(I) 

OP TO 90 

'88 on 89 I = 1",'M0 - V“-r- 

89 HI!) = -P. (I) + FLAST*HU), 

90 CONTINUE . 


v'J =1 .. 

9? ‘ A H ( JJ ) =' 0 ."0“’00 *’ 

on 95 KK a 1 , M J 

1 MS=I1(JJ,KK ) - 

JF(M8.F0 .0) OP TP 95 . 

AHfJJ)- •"= AH ( J.l ) +A (;JJ7:k(C^H(MST 

95 CPNT 1MIJF 

TJ = JJ + 1 ' 

I F ( JJ.GT.MF) c-n TO 97 JT 

GO TO 9? " - - 

97 CPMT ! hill F 

'CALL PATSliM (AH, 6'E iKMMf AHTAH ) 

JF( AHTAH. LF. 1.0 0-1? ) GO TO 2300 

' ' F L A M J = 0.0 00 •. “ + Tr}' : 

HTR a 0.0 00 

on 3). J a 1% Mi l. - 

31 HTR = HTR + HU)*R(J) 

ALPHA"’ = 0. O'DO" " 

FLA M j = -(HTR /AHTAH) . 

IF ( FI.MLST. c o'.o'. : o' hnF^'W’ : 96'' 

ALPHA = FLAHJ/FL'-’LST 

I’FIKnnNT.FO. 2) ‘ GO TO ~°3 

on TP 96 

' 93 HRTTF(6,617) KOI I A]T » ALPHA O'" 

96 FLFL8T = FLA’-.I '• , 

on mo i = i t mi.i • -'-••o — “ -- ■ 

ion x ( j ) = x ( ] ) +fla«j*h ( i 

' ' I F ( KOLIMT. GF ,W|I ) GO TO -106 ' ' “ *r~ ““ 

GO TO 105 ' 

106 00 104 'I = 1,MU,3~ T“ ; • ~ 

JF(DA8S(XLST( 1 + 1 )-X( 1+1 ) )/. LF . 1-. 00-07. AMDvOAPS ( XL5T t 1+?)- 

l.LF. 1.0 n-07 ) on TO 104 • ™.- " - “ ' 

f-P TO 105 

104 CONST I •'•'I. IF ' •' - r -~- -~ 

HP1TF (6,606) KKOUMT . , 


HR JTF( 6,6)3) 
WP1TF (7,61? ) 


R T R, K Oil m T " 
OTP , KOliMt.p;' ; 
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WR I TF ( 6 , 6 1 7 ) . K 0(1 MT , A L PHA : 


X( I + ?) ) 



vpitf ( 6,6 on* ) mk 

MR I TF f 7, 7006) NK, K PUMT' ■ 

US =1 

no 30? I = 1 , .•••! 1 , 3 ■ ■••• • ■ 

MM TP (6, 6600) JS,X( I )_, X (M+J. ) , X U +? ) , K° 

?o? cn T JMHF 

I F ( KK.OUMT . FO. 5 ) fen ' TfT * 

GP TP 19 

18 "R TTF{ 7,613 i RTR,KfH!MT " * . 

MM TP (6,6006) MK. 

WR T TF ( 7 * 7006 ) MK'rKOUMT ~f” 

*'5 = 3 

DO 303 " I "= ' 1 , Ml) ,'3 ••• 

IrfRITF (7, 7600 ) JS , X ( I ) , X ( 1 + 1 ) , X (I 4 ? ) , KO 

JS = JS + 1 

303 C.PKT 1M.IF 

19 KKOUMT = KK DIJMT + ~ — 

) F ( K KOU MT • PE • 6 ) GP TO 17 

GO TO 105 " ;■ ■ 

17 HPITP(6,606) KKOUMT 

GO TO 1 ; 

105 CPMT1NIE 

on - 107 i = l , mi i — ' 

•107 R l. S T ( I ) = R ( J ) ■ 'I 

no no i = i,mij v 

no r m = o.o no 

10? on 101 K = 1 , M I 

m s = 1 1 ( j , k ) - - 

101 RIMS) = RIMS) + F L A M ,1 * A I J , K ) * A H ( J ) 
IFIJ-MP) 102, 10?, 10? 

103 CONTINUE! ' . — 

rr- lop i = i , mi ) 

108 R ( I ) ■= RII.) + RLST(I) 

CALL PATSI!M(R,MII,KM,RTP) 

A 8 C0NT I Ml IP ■: ••••: 

UR 1 TP ( 6, 60? ) RTR , KPLIMT- 

RTH = 0.0 00 ; ■ -r : ~~ 

DP 91 J = 1,MM . 

91 RTH = RTH + R(J)W(J)T~“ — 

8 A I F ( DA 8 S ( RTH ) . Gt . 1.0 P-09) GO TO ?01 

GO TP 8? - - — ; - 

?0). URITF (6, 60? ) RTHfKnUMT 
'8? comtimiif — f -~ 

RLSTR = 0.0 DO 

on u 7 ' i = i , mu — - — : - 

A 7 RLSTR = RLSTR + RLST < 1 )*P ( I ) . 

IF (PARS (RLSTR l.GT.l.O n-09 ) ‘ GO TP ?00 

GP Tn A9 

200 HR 1 T F ( 6, 60 A ) ' RLSTR , MOUNT ; - :•• -•-•• t - •• • ~ 

A 9 rip 50 1=1 , Ml! 

50 XLSTC I ) = X I T ) — 

JJ =1 

' 35 AR( JJ) = 0.0 no — 

DP ?3 kk = !,MJ . 

MS = 1 1 ( JJ , KK ) . - ----- 

IF(MS.FO.O) Gn Tn ?3 

~ AR(JJ) = A P ( J J ) + ' A ( J j; KK'IPP (MS) ' 

33 CPr T } M IF -bu ~ 

'JJ = JJ + 1 — “••: 



>o 


■ )F! Jj;.0:-T - : <■■ 

" on to 35 

3 5 G POT T Mi.! F . ■/■& v 3 ' V ./ ■' . 

"' ARTAH •= n .O'nfj-'.'”' 
pn 36 1 = lff'T- ■ ■ 

36' 'ART AH' ="■ ART AH"’V ARTil CAH'CI 1 "T"' 

J F { KPi'M.T .FP.U on fh" 115 

113 I F (,RTR . LF . FN'i IU ) GO TO' 2)00 7'. \ 

J F ( KP’UUT - ( K DU’ 'T /.!•).)* U ) 114,115, 11/. 

1 1 5 HR I TF ( 6 , 6'onfc 1 KOIIMT ’ "" '' 75 * •" "tv? 
HR 1TF,< 6, 60?) P.TP,Kf.UMT ‘ 

W p. ITF ( 6 , 60/- T. ) ' R t F T R" V’’’ 

V! R I T F < 6 , 60 ? 1 ) R T H ' V 3 : ." ; ’ 1 7 • 

W R I T F ( 6,6400) Kr>UHTVAHTAH : t"* M 
V R. 1 T F ( 6 , 6 3 7 ) K DO UT , ALP HA , 1 

mk’ :r tv - 

WK , K HUNT 


,mu,3 . ■: •• . ■■ , r ■ 

j s’, x ( i ) , x { i +i i , x n+2) 


30). 
1 1 4 


100 

300 

2000 


’’WRITF! 6,6006) 

HR I 7 P ( 7 , 7006 ) 

”J F '~T 

P 0 301 T = ) 

W P I T F ( 6, 6600) 

OF = JF + 1 v :3 'O-v.x-3 •. • .... 

COMT li'iDF ’ - 

CPMTlM'F . “ '•••'• 

■ k niiMT = k nil mt • ■ + ;•■ r ~^yr 
] F ( K PI ifT . F.o. IMF 6 ) . OP TP ?000 

if(kpupt.i;t .hfa)X -on' jn 40 “ 

V P I T F ( 6 , 611 ) . KOUNT , R T RO 1 V- 
WP.TTF { 6,614) K 01 j UT;fA PIT A H 
WRIT F ( 6,6 10 ) KnUAiT HP.TR. r . : 

HR I TF ( 6, 6041 ); V’R LF 
HR] TF( 6, 603 D OTH: ' 1>57- 7. ’ •- 1 

HR ITF ( 6 , 6400 ) ’ KOUM.T , r AHT AH 
KOUMTVp l PHA 

MK ,KOUUT 


300 

2 600 
602 
603 
'604' 

60 6 

610 


WPITF ( 6* 61 7 )\ 
"HR ITF (6, 6006) 
HR ) TF (7, 7006 ) 

jf =' i "". ?::■ 

DP 300 I . = , 1 
UR I TF ( 6 ,6600 ) 
WRITE (7,7600 ) 

'JF': = ''.IS + "i 

COUTH UP 
"on to i 
URJTF (.6, 6100.) 


, MU , 3 

JF v X ( I )Vx ! I+lt, X ( I +7 ) , KOr 

J F , X () ) ,’X ( I + 1 ) , X ( I +?. ) , KD ' 


l.( J ) , LKj’J ) ,M(:J) , MK ( J ) , 1 1 ( J J .V’KI. ( J ) , J 


FORMAT (1H , • RTR= ' ,025. 16', 1 OX, ' KOUUT=- .' {I S K ’ 

FORMAT! 1H , ' R A n CO.U Jt.'OAr.Y RTH ,D25 . 16 i '■ KPUM 

FORMAT ! 1H , 'RUST NOTIgOM “JUGATF Tol'Rl ,1'OX , 'RL’FTR = ' , P2 
l*KnuUT= ' , IF ) -s , , ;J ’ . _ - 

FORMAT (1H1, •■*#«**•■ “ SPixlTT I OM l/FCTOR- DOFS >fhf I-MOR OVF "*'* 
1 1? , 5X , ' T If FF ' ) - I/-- ; ; -1 ; V 

FORMAT ( IH1 , ' IT = R AT IOfiF COMPlFTFO ' , 20X , ',KOUMT= ' , I F, 5X, 


T = ' , 15) 
5.16, )0X 

•P.TP= 


) P25 . ). 6 ) ...o'- ■' : ' H:'- 

611 FORMAT (1.H1 4 , PTrH' VAr'.!I5H^;V*^P^T , ^ n, . , A'T= Th‘T510?5. 16) 

612 FORMA T ( lHl , ' COLUMN. VVF-GT.pR*^:'*! 15) ’:1" 

’613’ FORMAT ( IX , ' F P L H TT 0 FCTOR''. DnE S MOT. 'IMPROVE * , 5X , * P TP,= » , 01 T . 8 , 
1 2 X , •KPUMTa 1 , 14)1;.- .:;.' .’• o'.j; . . ■ 7; • 

616 FORMAT !).H1 , * T 5 R M-I;W A TF • iF IT R ' “ATTr'A H = 0 • '/ / * KOUMT- • , 15, lOXJ « AHTAH= ”» 

1025.16) . o; ImM % 

6 1 7 FDR M A T n H , * K Hi i MT = -’ *1 ’1^1 1 OX VP' A L Pw A =' OR 5V. 16 , • < < « n op = r ROUMO’ 

lCPMOTTinu mum I- FP. . «<<<<•, ) .= "j 

5200" FORMAT! 15) ' "" 1 7 ’ 

6000 FORMAT! IX, 2 14, F7.ilVlF^l5 , I3'V7 1431? ) 

6001 FORMAT ( 1H0, »: MU »:>/3l6 ) ’ : J 7" ; 

■7‘ -V V ■- ' ' 

- • ..-61- 7 . ’, 



A no (t P03MAT ( 1HO, ?X , TOVAR I ANCF VECTOR FOR ; COLUMN*’ ' * 15//) 

6031 FORMAT! 1H t ' RTH= ' , n?? i 1 6 ) ” ” ~ 7 " ~ ’ 

AO 41 FORMAT ( 1 H . , *pl<:TR= • ,025. 16.» 

MOO FORMAT (l H , 1 SO” FTB 1 NP UROMO • //71 5’) ” 

A/,nr, FORMAT (1H , 'KiniM.T* ' , I A , 1 OX * ' AHT AH = 025* 16.1 

A 6 00 ’ FORM AT ( 1H , 1 r > , 3 F ? 3 . . 1 6 4X , 12) ’ 

4900 FORMAT ( 1 HI f 1 ITERATION « , I? ) 

7006 FORMAT { l OX , ' COVARIANCE ’ VECTOR ’FOR’ COLUMN* ’«/, I 5 , ?X . « KPIJMT* ' , 1 ? ) 
7600 FORMAT ( I ? » 3 F?3 . 1 6 » 4X » I ?.) 

e-poo return - 

FNO 
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) , L ( 990 ) , 
(990), 

600 ) , 


f. MAIN PROGRAM 'OF CG-ME.THOO FOR SOLUTION VECTOR 

IMPLICIT R E A L R 8( A - H , 6-7 ) * '7 ' " - V~.~ 

n I MENS ION. A ( 990,5 ) , II (990, 5 ) , X (600) , FL1990) , 9(990) , AX ( 990 

1 M ( 990 )",LK (990) ,MK( 990 ) , L'N( 990 ) ,ML ( 990 ) VMM ( 990 ) , TL (99 0)', XL 

?LM< 991) ) , R ( 600 ) ,H ( 600 ) , AH (990 ) , VC HP CK.( 99 0 ) , P. T ( 1, 600) ,RLST( 

r y L ^ T ( 600 ) , R 1 ( 600 ) ,XHV‘( °90 )7AR ( 990 ) "7 - V " 

VARIABLE : DIMENSIONS IN THE MAIN PROGRAM 

I) "VARIABLE : MATA CARD FOR '.NO' ME N 1 r 

C Nil = NUMBER OF UNKNOWNS 

C. '" ME a NUMBER nr- Of’SERVAT JOM ECUAT joWS ' • 

f. NT = NUMBER OF COLUMNS IN INDEX MATRIX Jill 

RFAD( 5,5000 ) MJ,NF,N1" " ” 

5000 FORMAT (3 15) 

CALL SOL N ( A , 11 , X , FI. , V , A X , L , M , LK , MK , L N, ML , MN , I L , K L , L M , R , H , 

1 A H , VC HECK , RT , RL ST , XLST , R 1 , CHV , AR , MU , ME , M I ) 

stop : 

END 


C 

C' 


SUBROUTINE SOl.N( A, 1 1 i X , FL , V, AX , L , M , LK , MK , LN, ML ♦ MN 
1 A H , V C. H E C K , R T , R L S T , X L S T , R 1 , CHV , A R , NU ,.N E , M I ) 
SOLUTIOM OF OBSERVATION FpUATIONS BY CG-FVE.THOO 


J L , K L , L M , R , H , 


IMPLICIT RFAL*R(A-H,0-7) 

" D I ME N SION A ( MF , MI ). ,“ I T (N E , NT ) , X ( NU > ,'FLTME T, V( NF ) ", A X ( MET," L ( 
IM(NE ) , LK (NF ) ,MK(MF ) ,LN(NF ) , ML (ME ) , MN ( ME ) , IL (.NF ) , KL< -MF ) . L 0 
PR ( MU ) ,H( Ml.)), 'AH (ME ) ,V'CHFf.K (ME ) , PT ( 1% MU) VR'LSTI Ml! ) , XLST ( Ml) ) , 
3,CHV(MF ) , AR (ME ) . \ 

"REWIND "1 : - 

VARIABLES IN SUBROUTINE ARE : A' , IM ♦> I R , KM. . 

'' Dl'FF" = T .0 "D-D 7 r - ------- -• - 


ME ) ," 

(ME), 

R1 ( M(l) 


NU2 = 

2 9 MU 

N'E2'"‘s 

PR ME ‘ 

NF3 = 

3 R N F 

NE'4""= 

' 4RME 

ME 6 = 

6RNF 

~ FM'I N 

- 1.0 D— 0 9"* 


C 

C 

C 

c 

c 

c 


N = INTERVAL FOR WRITE-OUTPUT 
M = a 00 : 

IM = INTERVAL FOR CALCULATING REE 

in = poo 

IP = MUMBFR OF- CONST RAINS 

IR =-'13 ' ' " 7 - . 

IK HUNT = l ••• / . 

KKOUNT a J '• 7.'7f77 ; 


(STATEMENT MR . 11 5 ) 
(STATEMENT MR. 7?) 


KO = 
KP -a~ 
I REE 
TREE 
I RTR 
TRTR 


CODE 

. p 

= CODE 
.’ = BO 
= CODE 
-= 81 


FOR 


ADJUSTED 


station COORD.; ..PUNCHED OUTPUT 


F 0° • RFE- rONCHF 0 OUT PUT 


FOR RTR-PUNCHED OUTPUT 


■ KM 'V PATSUM BA S I C"BTPCKS T-ZE"' ( FOR RTR );' — 

KM = 50 

F.KTR • =• 0 FLOAT ( ME ) /DFLOATTW I " ' r’—~ - 

K MM r: I D I NT < FKTP 9np LOAT ( KM ) ) 

INPUT : A-MATP IX , 1 1-NAT P I X "A NIV'L- VECTOR 

00 2 J a l,Mg . ■ v. 

READ ( 1 ) "7 L( J.).,M (•»?)♦ ( A ( J, K ) ,K=V, C ), I L ( J) , LK( . 

1 MK ( J ) , ( I 1 ( J , K ) , K = 1 , 5 ) , K L ( J > , M L ( J ) » MN ( J ) , FL ( J ) , L ' ( J ) 

? CPMT IMUE -*• 

DO JO J = l,vF . . v .. 

I F ( L ( J ) ; ME . LK ( J ) IV OR". I-Vrj ) .'NF .' MK CJ I .PR . I L ( JIT NF .KIT,')) 

. - 63 - . 


Go TO" 2500 





I F ( L ( J ) • M.F • F'L ( J ) . Oft . ( J ) . MF. MN ( j )•) r,fl TO 2600 

) 0 CF'NT } M.lr -cy.rr'-: ~ " 

WRIT r ( 6,6000 ) <L< J ) ,w (.)>, ( A( Ji. ; K) jX=l » c ) , T L ( J )', L K ( J > , 

] I'K ( J 1,11! ( J » K ) iK = l,5) , KL <J V,NL'( v ) , .’••N t ) , FL ( J Y? D'( J > , 0 = 1, WF'1 
WR I T F ( 6 , 600 ] ) MIJ t MIF»Ml . 'V55Zfi ;'Z / ' . 

CALL PSr.A.LF't l . o ■ n-0 2 » A V r : F i K.' T^.T. ~~ 

CALI. f‘Sr.ALP Cl.0 0-02., Fl 
non ;i =• i,Mu 

11 X(J) = o.n no • - .... • 

A FT PR THIS CARP READ X-VFCTOP' I F. MORE ) T.FR AT t OOS AR F WFFDFD WITH' GTVEM X 

12 contimuf v .. ; ; i .' ; 

Tie i (■ ; "**■*■-- ..... ....... 

16 R( J ) =. 0.0 00 __ ':■■■■ 

re* 15 j = 1 , n c- ' '• . "... • 

vr.HFCKU) = o.o no 

15 v ( o ) = 0.0 no ■ " ~ - - 

J - i . V'-7 . ' • 

"20 a y t .! ) = o.o no X“ 

DO 30 K. = 1,01 \. 

PS = I 1 ( 0 , K 1 - -- - ------ 

IF (OS. FO.O) 00 TO 30 ' : w- ; 

A X ( J ) = A X ( J ) + AQiK)«X|MS.)-#..f 
30 CONTINUE .. •vr'Z ., v 

j = o+i 

IF( J.GT.iVF) GO TO 40 . 

or TO 20. i^-rrr-fr-frr 

40 COMT 101 IE \a- : 1 

' CALL DOHA 00 ( A X »’ FL ,V', (5 .*’1' 1 '"j’T. 

VPV = 0.0, 0.0 

CALL P A T S U l-‘ ( V V W F , X pp. , VPVA) " 

CALL P'SC.A’L-F ( 1 .0' 0 04»VPV:,1»1)>j- - 

WRITF ( 6, 6070 ) .VPV 

F LA. ST = 0.0 DO 

FLKLST = 0.0 DO 

RFF = 0.0 00 

" Km i NT =' i z *~ 

RTR = 0.0 no 

45 RTRLST = RTR r 

DO 46 K =. 1 * Mil . 

XLST ( IO = X ( K ) "" ' ~ “ 

46 R L S T ( K ) = P (K) 

bo 50 K '=1',M(I ' 

50 R ( K ) = 0.0 DO 

= 1 ‘ 

65 DO 60 K = 1,NI 

"'"OS = 111 J»K) “■ ■' * 

60 R ( (•! S ) = R(f'S) + A ( J , K ) *'/ ( J } 

J = " J + l " " 

IF(J-MF) 65,65,70 

70 CTMT JHIF ' -- - .- 

I F ( IX Ot IHT— ( IKO| imT / 10 ) & 1 0 ) 71,72,71 

72 or 7? k = i + mi i; .-.j. 

73 R 1 ( K ) = O.o no 

74 nn 75 K= 1,MI 

PS = 11(J,KT ■■ 

76 R 1 ( *•’ S ) = R 1 (MS). + A ( .! , K ) *VCHFCK( J ) 

lF(J-NF) 76, 74 ,76 

76 CONTI HIF 

RFE = n.o no 
DO 7 7 I = 1,M(! 
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77 RFF = RFF + ( R < 1 )-R ] ( I ) )*.*?. 

' WR IT F ( 6, 6P : 1 i 1 knijMT » RFF t KOWMT ; 

RFF = 3*RFF 

’ “IF ( RT R . L F . P.F F ) GO' TP )4'" 

GO TO 71 

• ' — WR I T F ( 6 , 607 I KOUMT ,'R F F , R T R 

GO TO 9000 

71 CALL PAT SI If-'- ( R * Mil* KM'/RTP ) 

WR I TF { 6 , 60 2 ) RTRiKOL'WT, IKPLiMT 

1 f < i knijKi t— t jkoiiMT/nka 1 > a? , ar ,4?. -.;~ 

43 IF ( IK01INT- ( JKPIIMT/IK )*IM ) /.8 r 

- "“ 44" WRIT r- ( 7/701 ) I k'dOHTVRTP t’K n« 'MT', I RTR"’ ' — 

WR I TF ( 7 * 701 ) IKOUMT , RFF, KOUMT, IRFF 

48 RLSTR = 0.0 no / ~ 

DO 47 I = 1 , (VII 

47 RLSTR = RLSTR + RL'ST I I )*R( I ) 

IF(OARS(RLSTR ) . GT . 1 i 0 P-09) GO TO 200 

go to 49 - - — -r: . — — 

200 WRITE! 6,604) RLSTR VlKOUMT 

' 40 IFLKOUMT.FO. 1 ) GO‘UP'83' 

RTH = 0.0 no 

O'O' 91 '' J = TVMUT -A.-—- - 

91 RTH = RTH + R ( J ) (-J >_ _ 

83 RTH = 0.0 00 

P . 4 I F ( DA PS ( RTH ) ,GT . 1 . O"’ 0-09 ) GO TO 201'“' *.' 

GO TO 82 

201" WRITE (6, 603 ) RTH, j KOUMT ‘ • "■ 

82 IF(RTRLST.FO.O.O) GO TO 81 

iq ' = " i ' • " ' : : ~~ ' ~ 

35 AR(vU) = 0 . 0 00 

do ??" kk r r, Mr~ ~ — 

MS= 1 1 ( JJ.KK ) 

IF(MS.EO.O) GO TO' 33 P"""" - -- - 

AR(WJ) = AR ( JJ) + A< JJ,KK )*R(MS) 

33 CPK’T I HUE u *r 

JJ = JJ+1 

I F ( JJ.GT.f- F ) 00 TO ?.4 — 

GO TO 35 

34 CPN'T ] MU F ' ; ■' : 

artah = o.o no 

' DP 36 ' ~r = 1 *ME “ — " r r ~' 

36 ARTAH = ARTAH + AR(7)*AH(!) 

FLAST e “ART AH/'A'Ht Aft~.‘ .~~r • 

81 IF ( K01IMT-1 > 85,85,88 

'" '85 DP "86 I = 1 , Wi 1 ;■- --•• - rr ‘ 

86 H ( I ) = -R (I I 

88 00 89 I = 1, All 

-89 H(I) = -RU) + FLAST4HU) " ~ ; — 

90 CONTINUE 

C OP T DOOM A L ITY TFST BET . RL'ST AMO, P AMP’ RV/AWO H — ~ 

HTH - 0.0 n 0 

• DP" 87 I ' = , Ml) — -a- A ; 

87 HTH = HTH * HI I )’*H( I ) 

EMRMP ] = PSORT(RTR)" — 

EMRMR2 = DSortIRTRLST) 

■ - EMRMH = PS OPT (HTH) “ • ~ - 

0RTH1 = ( R L S T R ) / ( F NRMR 1 ) * ( E MR MR 2 ) 

— — PRTH2 -= ( PTH ) / ( FMRMR] )'M FMRMH ) 

IFfDARSIORTHl J.GT.liO 0-05. OR .P AR S ( OR TH2 j , GT . 1,0 0-05) GO TO 98 

- - - - - 65 - ^ 



°P HKITFI61WOI ORTH1 ,0RTH2 , IKOUNT 
99 JJ - 1 ' 

92 AH(JJ) = 0.0 00 

OP op kk = I , I ’ 

MS=Il(Jj.KK) 

i r- ( os . fo .0 j ; on to 9 r> 

AH(JJ) = AH( J.l)+A{ JJ,KK )*H{MS) 

'95 COOT I 01 1 E 

JJ = JJ + 1 ' 

1 F ( J J . c- T . 0 E ) GO TO” 97“ ~ 7 ‘ 

GO TO 92 

97 C, POT I 01 ip'”’ r .- ; - 

CALL PATSIIf-M AH, MF. , Kl'M , AHt AH ) 

1 F ( AHT AH.LE.l .0 0-1? ) "GO T0’'2?00 ' 

Fl.AMJ = o.n no 

htr ="o. n no ■ - ~ 

no j = - 1,01.1 

'?•]" HTR = HTR' + H( , |)*p.-( vll " 

ALPHA = 0.0 nn 

FLA M J = -(HTR/ A UTAH) ' --- - 

)F< FLMLST.FO.n.n on) Gn to 96 

ALPHA = FLA f-'J/F L M L ST "V " 

51 IF ( KOIIMT. F0.2 ) GO TO 93 

1 F ( IK HINT- ( I KOI.) NT/. 50 >*50 ) 94,93,94 

93 K'RITF-I A, 617) T K OIJMT * ALPHA • 

V P. I T F ( 6 , 6 2 2 ) KOHMT , FLA ST, FLAM J " ' 

94 COOT TNIJF 

C- 0 TO 9 A ’ " ; " • 

96 FLOLST = FLAM J 

no mo i = i.Mii" - -r 

100 X ( I ) = X ( I )+F|.AMJ*H< I ) 

I F ( JKPIIMT.GF.NII2) GO TP 107" " r 

GO TO 103 

107 no m2 i = i,‘MD ,3 : 7 ‘ 

JF( DABS ( XLSTf 1+1 ) -X ( I + 1 ) ).LE. 1 .00-04 . AMP. DABS ( XLST ( 1 + 2 ) -X ( 1+2) ) 

l.LE.1.0 n-04) GO TO 102 '"O' ' • - 

GO TO 103 

102 Of NT Jr UP " ' ■ T ~ 

WRITE (6, 606) KKCUIMT 

W P I T E ( 6 , 6 1 2 ) J k 01 INT , RTR - - — ; ~ “ 

WRITE? 7,701 ) JKOIIMT, RTR.KOIIMT, I RTR 

W P I T F. ( 6 , 6 0 J, ) -J KOIJMT , R FF *KOVNT\ - • ™*~ 

WRITFf 7,70). ) JKOIIMT, P.FF , KOHMT , IRFF 

WP.JTF (7,712) I KOI. I JIT , RTR " 

WRITFf 6, 60M ) RI.STR 

VRITF (6,6031 )RTH ‘ ~ ' 

WRiTF( 6,622) IKOL!MT,FLAST,FLAM.J 

VRITF (6, 6400) JKOIIMT, AHTAH " 

WRITE (6, 617) IKOHMT, ALPHA 
WRITE (6, 6006) 

JS = 1 

nr ?o: j i = i , mi i , ? • ~r — ~ 

WRITE! 6, 6600) . ! S , X ( 1 ) , X ( I ■+ 1 ) , X ( 

VRITF (7, 7600 ) JS,X('I ) ,X( I + ll;,X ( r+2) ,:<n 
JS = JS + 1 ' 

303 CONTINUE ! — — - — 

no 304 i = !,nf ; 

306 V C I ) ^ A X ( I ) + FL CD"” , . 

WRITE (6, 6007) ' 

WP ) T F ( 6 , 67 00 ) ( V ( I) , I = r, r ,: F -.7 — -y"" - ■ 

CALL P ATS I IM ( V , me , KMM , VPV ) c ; ■’ 

CALL (SCALE ( 1 . n 0 04 , VPV, 1 , 1 ) _gjg_ ' "f" 



WRITE) 6, 6070) VPV 
’ FM = ' DSPRT t VPV/ ( NE-NO+1 R ) ) 
WRITE) 6,60? ) FM 

KKPUNT = KKPI1MT "+ I “ 

IF ( KK0UMT.GE.6) GO TO 17 

GO TP ' 107 ~ J ~- -"r— 

17 WRITE) 6, 606) KKOUNT :v / 

GP TP '9000 •' 

103 CONTINUE 

IF) KPUNT. EP. 1) GP TP ' TOT ~ 

GO TO 104 

10 T PP 110 I = T ,ME " "~ 

1 10 VCHECK) I ) = V) I ) 


+ FLAN J* AH (I ) 


." GO TO 106 

104 DP 105 1= 1» ME 

105 VCHECK) IT = VCHECK { I ) EL AMJ4 AH ( T ) 

106 VCHCKV =0.0 DO 

CALL' "PAT SOM) VCHECK , NE : ,"KW", "VCHCKV " 

CALL MSCALEd.O D 04 , VCHCKV, _1 , 1 ) 

1 20 AX) JJ ) = 0.0 00 

DP 1?C KK‘ = 1 , N I ^ * 

MS = I1(JJ,KK) 

- 1 F( MS .EO.P) GP' TO" !? O' — 

AX(JJ) = AX(JJ) + A) JJ,KK )*X(MS) 

' 170" CONTINUE "" 

JJ = JJ + 1 . V ; 

I F ( J J . GT .ME) GP TO T40;T.‘ — 

GO TP 120 . 

140‘ CONT I Nil E ' 

VPV L S t = VPV i 

DP ' 150" I = 1 , N E, * ' - r~ ~ 

150 V) II = AX) I ) + FLU) ■ 

V P V * = 0.0 00 - •— — 

CALL PATS UM ( V , N F , K MM , V PV ) 

“ CALL N SCALE ( 1.0 0 04VVPV ,'1,1) 

CHNVPV = OARS ( VPV-VPVLST ) 

*1F(CHNVPV.LE.1.0 "D-0‘4').” GP fp 2200 " ' 

GO TO 119 . 

2 2 0 0" V R T T E ( 6 , 6 1 8 ) K PUNT, CHNVPV " : " ~ 

119 CONTINUF 

DP 116 I = 1 t ME 

CHV(I) = DA PS ( V ( I ) - '/CHECK (ID 
J F ( CH V ( I ) .GE . 1 .0 D-02 1 “T-P TP 1 16 ‘ * 

GP TO 117 

116 "WRITE ( 6‘t 67.7) T,V) I) , VCHECK ( I ) ~ * ‘ 

IIP. CONTINUE ■ p 

- 117 CONTINUE — : — : 

I F ( OARS ( VC.HCK V-VPV ) . GT . 1 . 00-0? ) GO TP : 2 POO 

GO TP 11? " " 

2800 WRITE) 6,616) KPUMT , VCHCKV, VPV : : 

11? CONTINUE 

IF ( RTR .LE.FMTN) GO; Tfl ; 2T00 " V • 

IFIRTR.LE.V.O n-OR) T-P TO '115 

IF ( IKOUMT.FP. 1. ) GO TO 1.1 5 ' 

if) IKPUMT-) iKnuNT/M)*Nr ; ii : 4-,- ii5vn'4":~" 

115 WR I TF( 6,6000! KPUNT, IKOUNT 

* WRITE (6", 60?) RT.R, KPU NT VIXPUNT - 

WRITE) 7,701 ) IKOUNt, RT.P, KPUNT, IRTR 

" WR I TE ( 6", 601") 1 KPUNT ERE E'rKPUNT "" 1"-= 

WRITE) 7,701) IKPIINT, RFE, KPUNT, IRFE . 

" WR IT F ( 6 , 6041 ) r LS TR' r^r-^r vT • - v- - 

- 67 -. 



! HR ITE ( 6, 60? 1) RTH 

! WRITER, 622 ) I KOI i w T » Ft: A" $ T* FLA M J" " 

WR]TF{ 6,6400) I KOUNT- , AH'TAH . 

WR ] T F ( 6, 6 1 7 ) I KOUNT , ALPHA - T — ~ 

WRITE! 6,6006) ' 

on 301 i - i , mu , 3 •' 

' WRITE (6, 6600 ) JS,X( I ) ,X (TLl ) ,X (T+?) ,K0 
JS = JS +1 ! 

30] cr^Tiwf- " •'* : — 

WRITE (6,6007) 

wr iTF (6, 6700 > "( vn ) ,T=r,N'ETr”"~^ 

W R I T F ( 6,6070) VPV 

hr i t f ( 6 , 60 RO ) vr.HCK v “ ' : :~t 1 ~'-: 

FM = OSORT < VPV/ (NE-MU + IR LV- 

WP.ITF (6, "606) FM • — *--- ~~~ T 

116 CONTINUE 

I F ( I K OU NTt ( 1 K 01 IMT / I M ) *T ’’ T^' TQ , 80,7 9 
80 IF(RTR.LF.PFE) - GO TO 78 . " . 

78 WRITE! 6,613) IKOUMT , RTR , RFE, KOUNT 

I K 01.1 NT = 1K0I.INT + 1" ' “ 

GO TO 12 ... . 

79 KOUNT"= 'KOUNT’ +'T — 

I K OIJNT = Ik 01 1 NT.' +• 1 

1F( IK0Ut'T.FO.MF6) GO TO "2000 — 

IE (KOUNT. LT..NE6) GO TO 65 

2000 WRITE (6, 61 OV 1K01.INT' ~ — r- 

2300 WRITE! 6,616) I KOUNT , AHT AH 

2100 WR’ITF (6,611)’ IK Oil N T , R t R ~ ' : '~T 

WRITE! 7,701 ). 1 KODMT, RTR , KOUNT. , I RTR * 

HR I T F ( 6 , 60 U I KOI ll-'T , P FF , KOUNT " 

WRITE! 7,701 ) IKOUNT ,RFF, KOUNT, IRFE 

WRIT E ( 6 , 6061 ) RLSTR — - : -- 

WRITE! 6,6031 )RTH w- ! 

WRITE! 6* 622) I KOUNT , EL AST* FLAM J ’ 

WRITE! 6,6600) I KOUNT ,AHTAH ' • ' 

WRITE ( 6, 617 ) 1 K OU NT , AC- P W “ ™"“ T 

WRITE!. 6, 6006) . - v - ^ " .L 

no 300 1 = 1 ,NU,3 

HR I T F ( 6 , 6600 I" J S , X < I ) ,'X'(T^n *X ( I.+ 2 ) , K 0 

WRITE! 7,7600) .1 S , X ( I ) , X (I +1) fX ( I + 2 ) , KD_ 

300 CONT IMIIE 

WRI T E ( 6, 600? ) ' ~ 

CALL jiSCALE ( 1 .0 0 02 , V , ME ,"). ) " 

-WRITE ( 6,’ 6 700 i ' (VI 1 1 ',T=f, MPcF--~"T ~~ — — ~ 

WR I TF ( 6 j 60 70 ) vow 

WR I T F f 6 , 6080 ) VCHOKV ;■ -y^r ■ . 

FM = DSORT ( VPV/ ( NF-NU + IR ).)',' L , 

GO TO 9000 - 

25.00' W P. I T E ( 6 ,6 100) L ( J ) VLfCt J TfM (7J )“f'N K { J ) V I L ( J ) t K L ( J ) , J— 

2600 WRITE! 6,6.200) L ( J ) , ML ( J ) , M ( J ) , HN ( J ) , J 

"601 FORMAT! IF , I KOIJM.T = ,1 » R PE ; = ”• , OPE'.'I 6 , 5X K. Pl.t NT's ' j'5T”“ 

602 FORMAT ( 1 H , ' P TR = • , 02.6^ |.'6 ,-10X , • KOUNT= * , J? , 10X , * I KOUNT = » , 151 

60? FORMAT! 1H , • PAO .CnNjtTGA-.nM^f'.'RfB*' V n ?5. •Tfe-,'' 1 IKO!INT= ' , I 5 ) " 

606 FORMAT MK , 1 R L A S T NOT: CM>$VGATF -TP R ' , I OX, 'RLSTR = » , 025. 16, 10X, 

605 FORMAT! 1 H , 'S T A MOA RO : ,Er4dM OF UNIT WEIGHT.#, »,ni2.6) 

606 FORMAT! 1W1 , ■ * •-;:<< * 6 S OL-UT- 1 ON' V E C. TOP OOFS NOT IMPROVE **»**• , 



115, 5X, •TIMES' ) . , 

607 ' FORMAT ( 1 H , ' T F RM I NAT JON OF I TR R AfjONS~ OHF TO P TR . I F. RFE *7/ 

1<IK0UNT= ' , I5,5X : V , !rF.F = ; ',0?5. U<,5X., 'RTR= ',055.161 

'610 FORMAT! l'hl , ' ITERATION? T. OM P L F T F D ' .', 20X7 ' I KOUNT = ' ',15V " 

611 FORMAT (1 HI , ' RT.R.’, UAnTSHFS ' , POX , »IKO0NT= ',!5,5X, 1 R TR = 1 ,025.16) 

61 2” FORMA T ( 1 K* y S 0 L" LET -ION ‘ ’ V F C TO R OOF S' MOT' ' I MPROVF • , VOX , • I KOHNT= ' , 15," 
15X , • RT-R= ',025.16) '• 

61 3" FORMAT ( lH'T ' IK0llNT = ' ,'I5,5)C, 'RTR='TD25 ,"16 ,'5X , ' R C F= ' , 025. 16', 5X , 

1 • K OMMT - ',151 ' V 1 ' J ' 

61 4 FORMAT! 1.HI , ' TF R Ml N’A f.F~'-*FPfi~"'A FIT AH=0 ' / / • TkOlIMT ='* , 15", if)X,' ' A'HT AH= • 

1025.16) , 

' 616 'FORMAT! 1 H " , 1 VCHCKl/.'NF.'V.PV » V 1 OX IT, 2025 .16) ' 

617 FORMAT ( 1H , ' I KOl'.NT='A;f 1,5 , 1 OX , ' AL PHA= • , 025. 16, '««IJPPER BOUNO 

' TCONDJTIOM Mi.iMP'FR'i * 1 - -- 

6)B FORMAT ( 1H , '*«’**** 0:"I'Kni!NT= ',15, 10X , ' VPV INCREASE* ',025.16) 

" 620 'FORMAT! 1H' , ' ORTHOGONALITY ' T FST : ' ' , ' 0RTM1 = R LS TR - ' , 01 6. 9 ,' 5X , ' " 

1 • 0RTH2 = RTH= 1 ,01'6, : 9;>5X, ' IKOUNT= ',15) 

622 '.FORMAT! 1H , ' “ KOI !MT : = T5V5X , • EL AST = ~ < ,T-?0. 10,5V, ' FLA M'J = •", F20.10) 

62 3 FORMAT (1H ,'I= ' , I 5;,' '■ , 'V = ' , F20. 10, 10V , ' '/CHECK* ' ♦ F20 . 10 ) 

701 FORMAT! 15 ,025noVfMWFX'n2 )'' " ■ * 

712 FORMAT ( *5X , » IKOOMf ^'>•'•‘♦^5 » • RT R= ' , F20. 10 ) 

6 0 0 O' FORMAT! IX , 2'I^T.F 7 . 2',;6'F 12.5,13, 714,13, 2T4,F 1 1.4,1?) 

6001 FORMAT OHO,' Nlj NF N I • / /3 I 5 ) 

'600 6 FORMA T ( 1 HO , 2 X , ' I * il2XY' 1 «y 22 X ,*'• PHI ' , POX V • LA MR 0 A ' ) 

6007 FORMAT ( 1 HO, » (.M&*** ' ) 

6031 FORMAT! 1H , ' PTH=* ,025.16) 

60 A 1 FORMAT (1H , ' RUSTrV>%025 . 16) 

6070 FORMAT! )H ", • VRV =7^^025 ) 

6080 FORMAT (1H , • VCHECKVi*;&025 .16) 

''6) 00" FORMAT ( 1H ~, * S0MF.ThMC'^!P'0M0 1 //7 15 T " '7"~ 

6200 FORMAT ( 1H , 'MOT MA’Tf.H.INO ' //5 I 5 ) 

6400 FORMAT! 1 H V IKpiMTs^I^-V 10X , ' AHT AH = '"»7D25 16' 

6600 FORMAT ( 1H , I 5 , ’3F23 6X , I 1 ) 

6700 FORM AT ( 6 F 20 . 1 03) ;T!T^y — - 

6900 FORMAT ( 1H1 , • TTER'AttlON ',215) 

7600 FORMAT ( I 5 ", 3 F 2 371" 6T5 fctT -IT : — — 

9000 RFTURN , . - : :j •, 

END 



' St ’PROMT I NF PATSUMt R 
IMPLICIT R, E A L * R ( A - H , 0 - 11 

DIMFMS'I ON "R C NIJ V " ' ' 

RTR = 0.0 no 

— jv>p t-tij/km-' ~ 

ML = MU~ i NP*KM ) . . 

I F ( N C • L F • 0 ) ~ G0"T'0' "2FT“7;ff^r 

ii s..(km*np> + I. ' '■' "i ff'-i 

~ RTRi ='' 0.0 00 ' — -ff-T” 

nn pa j. = ii ,Mi ■ ■ 

'2A'"RTR T' = -RTRI + R'Qi*R"FO 

RTR = RTR + RTRI. ■ T • A 

25 C ON T I (.'ll F "‘7 

IF (KM. OT. MID GO TO 28 
" DO 21 ' I = l , MP • 

II = I * K M A . 

- \!J ‘ I I "- ' (KM- IT-- 

RTRI = 0.0 DO 

"DP 22 J = JJVH 

22 RTRI = RTRI + R<j)*R(J> 

'21 RTR = RTR + RTRI • — ~~ 

RETURN 

28 * WR J T F ( 6 >'61'I r "' — 

61 FORMATtSX, ' IN PAT SUM, KM GT 

R FT IJRN " — 

EMI) 



"'SUBRnUT IMF DUMMY" A i B ? C t P i T t S I GW. iSCALE ,»N 1< ~~ 

DOlfBLF. PREP. IS ION bV ]< , W 1< , CP 1 < , P*1 < , m< f S I ON, SC A L E 

-'ENTRY NWRITF2A, IA,JA,MAMF< 1 ffl 

DJMFNSIPN N‘AME~3< • 

*WR ITF.r 6'i 1 <M ARE riXV I A ,- JA .* T~ 

1 FORMAT 5<A4,27X, I5,EH ROWSEX, 15 ,6H COLDSO 

" "J7J # I A * J A ' - --r- 

PO ? 1 1 r*. 1 , I A , 

2 WR'ITFK6V3<I 1 ,ffAZJ<V;i« ; TT,'j"J* I'A< ~~t • 

3 FORMAT?/'^! R0WJ5,//XlX,lP5020.fl«. 1 

RETURN 7 

ENTRY MPIINCH^A, IA, JA,NAN< 

- 'j'jifix ia*ja-i </3?;r •“ " 

HR]TE*7,4<5<A?fH*3-2< , AS JI*3-K, AX] I*3< f NA'O, 1 1 , II«1, JJ< 

- 4 FO’RM'AT* 3D24 . 1 6 V A'4', T'4< ; """” 

RFTURN , ! 

ENTRY -MMUL’T - A i’ R"» N » M‘» L’V’C'C ; 

ENTRY GMPRDXA t R,C,N,M,l< 

1K4-M 

: do io k«t,l ~~ v : 

.IKfclKCM 

] P #= I R f 1 _____ ' . • . 

] Bt.‘ I K '■ _____ _____ _ , , ' ■ y- 

DP lo;itfi*M ,___ ; _____ 

IPRIRf.l' h 

in cst-r <$ r.n r < fa ?.j i <*pt r r'R< ■ r — 

RETURN 

- -“ENTRY ’ GT P RO 7,A * R"V 0 VN » MVL'<' ~ 

N1#N 

ENTRY "T 1 Ml II T* A7WV M * B » Mi:VL7C< ~ — -t— - - 

] FS!N-N1<60 » 8 1 60 ^ __ ’ • 

jk#-n 

■ “ DO 20 K U 1 VL " i “ 1 7 

lv.o‘0 

IK«IKCN - : - *•••• — 

PO 20 J#1.M _ ' 1 _ ___ 

IRtfIRfil 

c»lR<so.nno — : 

no 2 o 

JRIMRSI 

••• 20 • ' 'CS?IR<«C^IR-<f ASI3<*R|?rR<-- •: 

RETURN . 

' ENTRY T2MULTSA~MtM-;P>L~»-M1 ♦ €< — 

1 F*M-I- 1<?0, 9, 50 

DP 30 K -i 1 * L 

' — - 00 *30 • JS1, N ' ■ — 

IP*K-L . 7 _____ 

UIRJ-N . 

— -py I R‘<;7 0 • 000.7 

• DO 30- I/U »M 



30 

SO ' 

60 

70 

100 


poo 

300 

A 00 


500 


600 



1 p.;. 1 i Rf,i. •;* <■. . . .. 

cr; iR<^c‘i'lR<'6A«.i r<*r*%I((< ~ v- 

R F T U R (\‘ . . 

M. .vM] •• ; 

■'OR I T F” 6 ,’70<M ;m"' “ r “ T^-f 

FOR [•■ A T Y 3 0 H 1 * * *N OT 6 VALID DPFRAT 10N-- JA6 T5 1P«T5, 

* 16H. CALL IGMORED.^- ~;V r— - 

P. FT URN 

ENTRY i-iSCALFOROALF , A V IT, J T<""' " - 

00?'; IA*JA . 

nn mo JS'I", JJ ” : 

t S' J < *'• S C A L E * A N J < \j. " - -■ ' " 

FNTRY f-'AOD^A ,R, I A* JA ,C< 

' S I OM;'J 1 . 000' “ “ ~~T~ 

00 TO ?00 

ENTRY" HR UR TO' A', R VI A , JA','T<" ~ T ~ ~~ ~ ~ 

S Jf5Mf.f-l.0po 

J.J;; I A* JA ' ” ““ " 

00 ?0o 0«1 , JJ 

cm j<«Ai-: j<f.s ir,N*p.K'j'< : ; : 

RETURN 

"ENTRY N0IIPRA,TA ,'J'A,R< 

0001 A* j A 

DO 400 J {, 1, OJ 

p ?• j < a a y. j < 

ENTRY .NATUAVA , I A, JA ,P,CiT< : " '.f, ' ; 

) K - 1 A " 

nn 6 no xor , ja — ; — . 

JK?; JKf. IA 

nn son 1 , 1 a ' • ‘ , ' ■ • - „ 

J R f ‘ J k __ , • ; . . _ ___ ■ 

T Y 1 T < 0 . 0 PO • 

DO 500 101, IA - . — 

I MO IN FI A 

T /’ I T < T Y I T < F- P Y I N< * A Y I P < 

nn 6 oo j •? 1 , j a - . 

CSjR<i‘ 0 .npo 

• nn son iff!, 1 A'- 

IR6IRF1 ' . 

c a . IR<:.‘C^IR<f AKIR<*T5:r< ' T “" 

R FT I.IRN 

"ENTRY MAMATSA , ! A V JA , P , CVT< — - 

IRrO 

'DO 800 K ?; 1 V I A :~TT “ 

1KRK-IA 

nr 700 j • 1 1 , ja 

I M," J— J A 

" TTOITFl ' — 

. TS']T<?;n.0pn 

nn 7oo' 1 a 1 , ja •••• ; — • — 

. - 72 - .... 



If 'i I' -'F.JA 

I'BjilBF. I A 

700 T 5< 1 T < * T 2 1 T < f, P V T «<# ft 2 1 P <r 

"T"' on POO JV'-l ,IA 

] P:“- J- I ft 

“■‘IRvIRf.'l* 

C2 ] R<if 0 . OPO 

DO 800 I 1 , JA 

]R? I PC I A 

POO ■ • c 2' I ft < ='( c 2 I R < f. 6 2 1 0< * T 2 1 < “ 

RETURN 

END 




REFERENCES 


Adams, Oscar S. (1930). "The Bowie Method of Triangulation 

Adjustment as Applied to the First-Order Net in the Western 
Part of the United States, " U. S. Department of Commerce, 
Coast and Geodetic Survey , Special Publication No. 159 .. 

U.S. Government Printing Office, Washington, 

Ashkenazi, V. (1967), "Solution and Error Analysis of Large Geodetic 
Networks, " Survey Review, No. 147 , pp. 194-206. 

Ashkenazi, V. (1969). "Solution and Error Analysis of Large Geodetic 
Networks, " Survey Review, No. 151 , pp. 34-80. 

Badekas, John (1969). "Investigations Related to the Establishment 
of a World Geodetic System, " Reports of the Department of 
Geodetic Science No. 124 , The Ohio State University, Columbus. 

Beckman, F. S. (1960). "The Solution of Linear Equations by the 

Conjugate Gradient Method, " Mathematical Methods for Digital 
Computers, Vol. I, edited by Anthony Ralston and H. S. Wilf. 
John Wiley and Sons, New York. 

Bomford, G. (1965). Geodesy . Oxford University Press, London. 

Bouchard, H. and Moffitt, F. H. (1934). Surveying. International 
Textbook Co. , Scranton, Penn. 

Deker, Hermann (1967). "Die Anwendung der Photogrammetrie in 
der Satellitengeod'asie, " Deutsche Geodatische Kommission , 
Reihe C, Heft Nr. 111 . 

ESSA (1969). "Precise Traverse Chandler, Minnesota to Moses Lake, 
Washington," Environmental Science Services Administration 
Coast and Geodetic Survey, Rockville, Md. , May 12. 

Foreman, Jack (1970). "Spatial Traverse: Scale for Satellite 

Triangulation, " Paper presented at American Geophysical Union 
National Fall Meeting, San Francisco, December 7-10. 


- 74 - 





Fox, L. (1965). An Introduction to Numerical Linear Algebra . 

Oxford University Press, New York. 

Gaposchkin, E.M. and Lambeck, K. (1970). "1 969 Smithsonian 

Standard Earth (II), " Special Report No. 315 . Smithsonian 
Astrophysieal Observatory , Cambridge, Mass. 

Geonautics (1969). "Geodetic Sate ll ites Observation Station Directory, " 
Prepared for National Aeronautics and Space Administration. 
Geonautics, Inc. , Virginia, July. 

Gergen, John (1970). "The Analysis of a Short Segment of the U. S. 

Coast and Geodetic Survey High-Precision Transcontinental 
Traverse, " Master of Science Thesis, The Ohio State 
University, Columbus. 

Ginsburg, T. (1963). "The Conjugate Gradient Method, " Numerische 
Mathematik, 5 , pp. 191-200. Springer Verlag, Berlin. 

Gossett, F. R. (1950), "Manual of Geodetic Triangulation," U. S. 

Department of Commerce, Coast and Geodetic Survey, Special 
Publication No. 247. U. S. Government Printing Office, 
Washington. 

Gotthardt, Ernst (1968). Einfuhrung in die Ausgleichungsrechnung. 
Herbert Wichmann Verlag, Karlsruhe. 

Grossman, Walter (1$61). Grundzuge der Ausgleichungsrechnung . 
Springer-Verlag, Berlin. 

Grossmann, Walter (1964). Geodatische Rechnungen und Abbildungen 
in der Landesvermessung . Verlag Konrad Wittwer, Stuttgart. 

Helmert, F. R. 41380). Die mathematisehen und physikalischen Theorien 
der Koheren Geodasie, I, Tei'l. B. G. Teubner Verlag, Leipzig. 

Hestenes, M. R. & Stiefel, E. (1352). "Methods of Conjugate Gradients 
for Solving Linear Systems, " Journal of Research of the National 
Bureau of Standards, VoL 49, No. 6 , December, pp. 409-436 

Hilger, F. & Remmer, W. (1967). "Das Rechenprogramm AUGL zur 

netzweisen Ausgleichung einzelner oder mehrerer Triangulations 
Ordnung in einem Guss, " Deutsche Go o dUtische Kommission , 
Reihe B, Heft Nr. 139. 


- 75 - 




J o rdan/Egge rt/Knei ssl (1958). Handbuch Per Vermussungskunde, Bd. IV, 
T . Tell . J. B. Metzlersche Verlagsbuchhandlung, Stuttgart. 

rlordan/Kgge rt/Knei ssl (1959). Handbuch Pe r Vermessungskunde, Bd, IV . 
2. Teil. J. B. Metzlersche Verlagsbuchhandlung, Stuttgart. 


Korhonen, J. (1954). "Einige Untersuchungen tiber die Einwirkung der 

Abrundungsfehler bei Gross-Ausgleichungen. Neue-Ausgleichung 
des sudfinischen Dreieckskranzes, " Mathematisch-Naturwissen- 
schaftliche Sektion der Philosophischen-Fakult'&t der Universit'dt 
Helsinki, 26. May. 

Krakiwsky, E.J. (1968). "Sequential Least Squares Adjustment of Satellite 
Triangulation and Irilateration in Combination with Terrestrial 
Data, " Reports of the Department of Geodetic Science No, 114 . 

The Ohio State University, Columbus. 

Lauschli, P. (1959). "Iterative Losung und Fehlerabschatzung in der 

Ausgleichungsrechnung, " Zeitschrift fur Angewandte Mathematik 
und Physik, Vol. X. BirkHauser Verlag, Basel, pp. 245-280. 

Meade, B.K. (1967). "High-Precision Geodimeter Traverse Surveys in the 
United States, " Paper presented at the XIV General Assembly of 
IVGG, Lucerne. 

Meade, B.K. (1969a). "High-Precision Trans-Continental Traverse Surveys 
in the United States, " Paper presented to XI. Pan American 
Consultation on Cartography, Pan American Institute of Geography 
and History, Washington, D.C. 

Meade, B.K. (1969b). "Corrections for Refractive Index as Applied to 

Electro-Optical Distance Measurement, " Paper presented to the 
Symposium on Electromagnetic Distance Measurement and Atmospheric 
Refraction, International Association of Geodesy, Boulder, June. 

Meade, B.K. (1970). Private Communication. July. 

Mithchell, Hugh C. (1948). "Definitions of Terms used in Geodetic and 

other Surveys, " U. S. Department of Commerce, Coast and Geodetic 
Survey, Special Publication No. 242. U.S. Government Printing 
Office, Washington. 

Mueller, Ivan I. (1964). Introduction to Satellite Geodesy . Frederick Unger 
Publishing Co. , New York. 


- 76 - 



Mueller, Ivan I. (1969). Spherical and Practical Astronomy as Applied to 
Geodesy . Frederick Unger Publishing Co. , New York. 

Muller-Merbach, H. (1970). On Round-Off Errors in Linear Programming . 
Springer-Verlag, New York. 

Pellinen, L. P. (1970). "Expedient Means of Joint Processing of Ground and 
Cosmic Triangulation, " Bulletin of Optical Artificial Earth 
Satellite Tracking Stations- USSR. Joint Publications Research 
Service, Washington, D.C. 

Rainsford, H.F. (1955). "Long Geodesics on the Ellipsoid, " Bulletin 
G^od^sique, Na. 37, pp. 12-22. 

Ralston, A. (1965). A First Course in Numerical Analysis. McGraw-Hill 
Book Co. , Mew York. 

Rapp, R. H. (1969). "Geometric Geodesy Notes, " Unpublished. 

Schmid, Hellmut H. (1965). "Precision and Accuracy Considerations 
for the Execution of Geometric Satellite Triangulation. " U. S. 
Department of Commerce, Coast and Geodetic Survey, Rockville, 
Maryland. 

Schmid, II. H. and Schmid, E. (1965). "A Generalized Least Squares 
Solution for Hybrid Measuring Systems," U. S. Department of 
Commerce, Coast and Geodetic Survey, Rockville, Maryland. 

Schmid, Hellmut H. (1966). "The Status of Geometric Satellite Triangu 
lation at the Coast and Geodetic Survey." U. S. Department of 
Commerce, Coast and Geodetic Survey, Rockville, Maryland. 

Schmid, Hellmut H. (1969). "A New Generation of Data Reduction and 
Analysis Methods for the Worldwide Geometric Satellite 
Triangulation Progr'am, " Paper presented at the Department 
of Defense Geodetic, Cartographic and Target Materials 
Conference, October 30. 

Schmid, Hellmut H. (1970). "A World Survey Control System and its 

Implications for National Control Networks, " Paper presented at 
the Canadian Institute of Surveying, Halifax, April. 

Schmid, Hellmut H. (1971). Private Communication, June. 


- 77 - 



Schwarz, H. R. (1968). Numerik Symmetrischer Matrizen . B. G. Teubner, 
Stuttgart. 

Schwarz, H. R. (1970). "Die Methode der konjugierten Gradienten in der 

Ausgleichungsrechnung, " Zeitschrift fur Vermessungswesen, No. 4. 

Shedlikh, M. (1970). "Role of Satellite Geodesy in Further Development of 

Continental Astronomical Geodetic Networks, " Bulletin of Artificial 
Earth Satellite Tracking Stations-USSR . Joint Publications Research 
Service, Washington, D. C. 

Simmons, Lansing G. (1950). "How Accurate is First-Order Triangulation ?" 
The Journal, Coast and Geodetic Survey, No. 3, April, pp. 53-56. 

Sodano, E.M. (1958). "A Rigorous Non-Iterative Procedure for Rapid Inverse 
Solution of very long Geodesics, " Bulletin Geod^sique, No. 47/48 , 
pp. 13-25. 

Sodano, E.M, (1963). "General Non-Iterative Solution of the Inverse and 
Direct Geodetic Problems, " Paper presented to the XIII IVGG 
General Assembly, Berkley. 

Stiefel, E. (1952). "tiber einige Methoden der Relaxationsrechnung, " 
Zeitschrift fur angewandte Mathematik und Physik, Vol.III , 
pp. 1-33. Birkhauser Verlag, Basel. 

Thomas, G. B. (1960). Calculus and Analytical Geometry . Addison-Wesley 
Publishing Co. , Reading, Mass. 

Veis, George (i960). "Geodetic Uses of Artificial Satellites, " Smithsonian 

Contributions to Astrophysics, Vol, 3, No, 9 . Smithsonian Institution, 
Washington, D. C. 

Wolf, Helmut (1950). "Die strenge Ausgleichung grosser astronomisch- 
geodatischer Netze Mittels schrittweiser Annaherung, " 
Veroffentlichungen des Instituts fur Erdmessung, No, 7 , Bamberg. 

Wolf, Helmut (1968). Ausgleichungsrechnung nach der Methode der kleinsten 
Quadrate . Ferd. Dummlers Verlag, Bonn. 

Wolfrum, O. (1969). "Iterative Verfahren der Ausgleichung nach vermittelden 
Beobachtungen und einige Beispiele ihrer Anwendung bei geodatischen 
Lagenetzen, " Deutsche Geodatische Kommission, Reihe C, 

Heft Nr. 143. 


- 78 - 


Zurmuhl, R. (1958). Matrizen . Springer-Verlag, Berlin. 

Zurmuhl, R. (1964). Matrizen und ihre technischen Anwendungen, 
Springer-Verlag, Berlin. 


- 79 - 



